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ABSTRACT 


Fundamental results developed by Wiener in the 1950’s are combined with new 
work in the area of higher-order statistics to develop and explore a general model for 
nonlinear stochastic processes. The Wiener model is developed for discrete nonlin¬ 
ear systems and its orthogonality properties are analyzed to characterize its output 
statistics. An efficient structured procedure for computing the fc*^-order statistics of 
the model output is formulated in both the time and frequency domains. Explicit 
formulas that exploit the structure of the Wiener model are given for computing the 
cumulants and polyspectra. A necessary condition for a discrete random process to 
be representable by the Wiener model is discussed. A computationally efficient pro¬ 
cedure is given for matching the model output cumulants to estimated cumulants 
for a given process by minimizing the squared magnitude of the error. Examples of 


applying this procedure to given sets of data are presented. 
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I. INTRODUCTION 


In the field of signals and systems, as in other scientific areas, we often have 
measured data that we wish to explain. This data is considered to be explained 
when we are able to develop an appropriate mathematical model which can produce 
data “statistically similar” to the data that has been originally measured. Another 
criterion for “goodness of modeling” could be prediction. When the detailed variation 
of the data is sufficiently complicated, or when we want to to model an entire class of 
data with different detailed characteristics, a stochastic model is often appropriate. 
The work in this dissertation develops and explores a general stochastic model for 
data which combines fundamental results developed by Wiener in the 1950’s with 
new work in the area of higher order statistics. The latter work was brought about 
by a desire on the part of researchers to move beyond linear Gaussian models and has 
been bolstered by new developments in high speed computation and digital signal 
processing. 

Although a significant amount of work has already been done, there are still 
many areas where researchers continue to work on building models to characterize 
stochastic processes. It is appropriate to cite some of this work (as well as some of 
the limitations) here. 

The stationary linear Gaussian model has dominated work on time series anal¬ 
ysis for a large number of years. In this model a time series (or discrete time random 
process) is represented as the output of a linear time-invariant system driven by 
zero-mean white Gaussian noise. The two basic assumptions of this model, namely 
the Gaussianity of the input process and the linearity of the system, has enabled 
researchers to develop a large number of procedures to identify the system param- 
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eters such that the model output dependably represents the data. When the linear 
system is represented by a difference equation of finite order, the model is called a 
linear Gaussian ARMA (AutoRegressive Moving Average) model of the stochastic 
processes. 

Among engineers, some of the most recognizable contributions to linear Gaussian 
modeling are due to Wiener [1]. While Wiener concentrated primarily on continuous 
time processes because of the need to process signals with analog electronic devices, 
others, such as Wold [2] and Kolmogorov [3] had similar results for discrete time 
processes. Wiener showed how the second-order statistics of the modeled process 
can be used to identify the linear system to within a complex scale factor. He 
showed conditions for which the complex power spectral function of the process can 
be factored to determine the transfer function of the linear system. In this case, 
if the modeled process is assumed Gaussian, the model is called the innovations 
representation of the linear Gaussian process. 

The class of linear Gaussian ARM.4 models has considerable limitations in rep¬ 
resenting stochastic processes. Among these limitations are the following [4, 5]; 

1. Linear difference equations do not allow stable periodic solutions independent 
of the initial state. 

2. The models are not suited for modeling data exhibiting sudden bursts of very 
large amplitude at irregular time epochs. 

3. The models are not suitable for data exhibiting a strong asymmetry (about the 
mean). 

4. The models may not be best for data exhibiting strong cyclicity. Nonlinear 
approximation for the regression functions may be more appropriate in these 
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cases, while the regression functions for linear Gaussian models are inherently 
linear. 

5. Since linear Gaussian models are inherently time-reversible, these models are 
not suitable if the data exhibits time irreversibility. 

It can therefore be seen that the linear Gaussian model is incapable of represent¬ 
ing a significantly wide class of processes. Since only Gaussian processes are com¬ 
pletely characterized by their second-order statistics, the higher-order statistics of 
non-Gaussian processes provide information that can be used to build an adequately 
dependable model. Recently the problem of modeling non-Gaussian processes has 
been closely associated with the higher-order statistical analysis of random processes. 
If the modeled process is far from Gaussian a number of approaches have been fol¬ 
lowed in attempting to build a reliable model. Among these approaches are ’’local 
linear”. Multi Variate Adaptive Regression Splines (MARS) and neural nets. In one 
general approach the model consists of a white non-Gaussian process driving a linear 
system, while in another approach the model is comprised of a white Gaussian or 
non-Gaussian process driving a nonlinear system. In addition, there are a whole host 
a methods that focus on non-Gaussian and or nonlinear processes of one form or 
another but are specific in the types of processes that they seek to represent. It is 
safe to say that none of the “general” and none of the specific methods overcomes all 
of the limitations enumerated above. We can however review some of these methods 
briefly. 

A considerable amount of work has been done in building linear non-Gaussian 
models for stochastic processes, mainly due to the availability of both analytical and 
computational tools for higher-order statistics and the relative simplicity of results 
that arise from the linear model. In such models the process is assumed to be rep¬ 
resentable as the output of a linear system driven by a A’‘^-order stationary white 
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non-Gaussian process with known statistics. For such a process the mean is constant 
and the other moments up to fc*^-order exist and are functions of only the lag differ¬ 
ences. Cumulants and polyspectra of the modeled process are fairly easy to compute 
and can be used to identify the transfer function of the lineax system to within a gain 
and phase ambiguity. The methods that have been proposed to identify the linear 
system in this type of models include both parametric and non-parametric methods. 
In the parametric methods the system is characterized by its ARMA parameters and 
the transfer function is expressed in terms of these parameters. Then matching a 
set of the model output statistics to those of the modeled processs enables one to 
construct a set of equations solvable for the system parameters. In [6, 7] Tugnait 
uses both the second- and fourth-order statistics to recover the poles and zeros of the 
linear system transfer functions. Giannakis and Mendel use the output cumulants 
to determine the AR and the M.4 orders of the system [8]. Fonollosa and Vidal use 
a linear combination of the output cumulants to identify the linear system in such 
models [9]. Several others have also made significant contributions. 

Non-parametric methods for identifying the transfer function of the linear sys¬ 
tem are also used. Lii and Rosenblatt [10] provide a method to estimate the amplitude 
of the transfer function from the power spectral density function of the process and 
then show how to estimate the phase from a ^’‘'•-order polyspectrum. Nikias [11] 
shows how to estimate both maximum-phase and minimum-phase transfer functions 
by applying a convolution between the third-order cumulant of the modeled process 
and the same order complex cepstrum. 

Although linear models for non-Gaussian processes have a significant analytical 
support, they are incapable of representing a large class of non-Gaussian processes. 
In fact, recent analyses have shown that at least from a theoretical point of view, the 
set of linear non-Gaussian processes is vanishingly small [12]. From a practical point 
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of view however it is advisable before modeling a process to test for its linearity. 
Tong [4] presents linearity tests of stochastic processes in both time and frequency 
domains. Tekalp [13, 12] establishes the necessary and sufficient conditions for a 
non-Gaussian process to be representable by a linear non-Gaussian model. 

A class of linear and nonlinear stochastic processes can be specified by the dy¬ 
namics of a system the output of which represents the modeled process. A significant 
amount of work has been done to specify such dynamic systems. For nonlinear pro¬ 
cesses the system dynamic representation (e.g., the difference equation) has a kind of 
nonlinearity according to which different classes of the stochastic outputs can be dis¬ 
criminated. In his book [4], Tong presents some examples of threshold models based 
on piecewise linearity. He also presents different types of autoregressive models with 
nonlinear autoregression and other types of specific models. Recent work by Stevens 
[14], Lewis and Stevens [15] and Lewis , Ray and Stevens [16] has applied Friedman’s 
M.4RS methodology to give a practical way of applying threshold models. 

Another broad class of nonlinear random processes can be represented as the 
output of a nonlinear system driven by a white random process. The input-output 
relations of the nonlinear system are described in the Taylor series-like representation 
with memory known as the Volterra series [17. 18]. The model is completely known 
when the kernels of this representation are identified. Powers et al. present appli¬ 
cations of using a finite-order Volterra representation to identify nonlinear systems 
[19, 20]. In these applications both the input and the output sequences are accessi¬ 
ble; therefore different order cross-moments between the input and the output can 
be estimated and used to obtain the kernels of the representation. In this application 
the Volterra representation must be finite. Other similar work can be found in [21] 
and references therein. 

The Volterra series representation is limited in its application because a rather 
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severe condition is imposed on the magnitude of the input process to guarantee 
convergence of the series. Further the individual kernels cannot be identified sep¬ 
arately unless the representation is assumed finite. To overcome these limitations 
Wiener proposed a model of nonlinear systems which is associated with an orthog¬ 
onal functional representation of random processes. He called this the G-functional 
representation because the functionals are orthogonal when the input process is white 
Gaussian [22, 17]. In the literature the structure and analysis of this model is based 
upon second-order statistics of the model output. The model kernels are identified 
via different order cross-correlations between the input and the output. So far, no 
one has proposed a method for determining the model parameters from output mea¬ 
surements alone as is the case when we try to fit a stochastic model to measured data. 
However, the use of higher-order statistics of the output suggests the possibility of a 
new' approach where this could be done. 

In this dissertation w’e introduce the Wiener model for nonlinear systems as a 
genera] model for a class of nonlinear processes that are necessarily neither periodic 
nor bandlimited. W'e also develop the complete details of the representation in dis¬ 
crete time. Wfiener’s original development of nonlinear systems was for continuous 
time and there are a few significant differences. Following this we provide a complete 
statistical analysis of the model. In particular, we develop expressions for the higher- 
order statistics of the model output in the time and frequency domains as functions 
of the model parameters. Finally we use these developments in an algorithm to iden¬ 
tify the model parameters such that the model output statistics match those of a 
measured data sequence. 

The remainder of the thesis is organized as follows. In Chapter II we summa¬ 
rize the properties of higher-order statistics of discrete random processes that are 
useful to our application. New' methods for finding cumulants and polyspectra of 
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linear and nonlinear functionals of Gaussian random processes are also formulated 
in this chapter. In Chapter III we describe the Wiener model of discrete nonlinear 
systems. The orthogonality of the discrete G-functionals is used to specify the model 
structure and an analysis of this structure is performed. In Chapter IV we formulate 
the model output cumulants and polyspectra and develop an efficient procedure to 
compute them. This comprises the main body of the work. The Wiener model has 
unique structure that simplifies the computation and which has not heretofore been 
exploited. In Chapter V we consider the problem of modeling a discrete random 
process using the Wiener model. We present the test of linearity and provide a def¬ 
inition of the class of discrete random processes that can be represented using the 
nonlinear model. We then give two examples of nonlinear processes resulting from 
the discrete Wiener model that match estimated cumulants of some measured data. 
The nonlinear equations resulting from this modeling problem are solved using the 
Extended Kalman filter technique. Finally, in Chapter VI. we present conclusions 
and recommendations for future continuation of this research. 





II. HIGHER-ORDER MOMENTS, 
CUMULANTS AND POLYSPECTRA 


Higher-order statistics and polyspectra can be applied in a wide variety of prob¬ 
lems in signal processing and system theory. By “higher-order” we mean the statistics 
of orders more than the second. Among higher-order statistics applications are the 
problems where either non-Gaussianity, nonminimum phase, colored noise, or non¬ 
linearity are of substantial interest [23]. In the theory of nonlinear systems, a system 
output can be modeled by defining the input-output relation using a Taylor series-like 
representation. Thus, computing an output statistic of any order involves the com¬ 
putation of higher order statistics of the input. Moreover, the output of a nonlinear 
system is in general non-Gaussian even if the input is Gaussian. 

A. HIGHER-ORDER MOMENTS AND CUMULANTS 

This section begins with a general set of definitions for moments and cumulants 
and then proceeds to discuss properties of special interest in this thesis. 

1. Definitions of Moments and Cumulants 

The /r-dimensional vector x = [xi, X 2 . • • •, is used here to denote a 
collection of real-valued random variables Xj for which the joint moments up to a 
sufficient order exist. Let u = [uj. U 2 -• ■ • i denote a set of real variables (i.e., 
u G then the moment generating function is defined by 

= (2.1) 
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where S denotes statistical expectation. The order moment 


= £ {x\X2 • • • ijk} (2.2) 

is the coefficient of the term j^u\U 2 • • • u* in the Taylor series expansion of (^)- 

general the moment of order m = i/i + i/j H- Vvk, namely 

f • • • Xfc*} (2.3) 


is the coefficient of the term 


in the same expansion. 


Vx\v2\--‘Vk\ 


U 


1 


u 


2 



The cumulant generating function is defined by 


<?i^(u) = ln<^;^(u) 


(2.4) 


and the order cumulant denoted bv 


CW = cum{j-ia’2---3”fc} 


(2.5) 


is the coefficient of the term j^UxU 2 - • • Uk in the Taylor series expansion of <?i>x(u). 
Likewise 

cum{x^x^ •••Xfc*} (2.6) 

is the coefficient of the term 

« 1 I ^2 “At 

in the same expansion. The A'‘^-order cumulant can be expressed in terms of the 
moments up to the order by the moment-cumulant relation 

(_l)(9-i) kl « 


* ^ n /'(l)'/ 42 )l . . . Z,(g)l 11 


n 


(2.7) 
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where the summation is over all the possible values of q. Similarly the cumulant- 
moment relation is 



E 


— n 


( 2 . 8 ) 


In particular, if the exponents are such that i^i — 1/2 - * • • — — 1 tlien from (2.7) 

we have r f '1 

c<.‘>=n *4-^| n ^4 ^.9) 

q lte«i J l‘€«, J 

where ki, K 2 , • * • and k, are partitions of the set of integers {1,2, • • •, k} and q is the 
number of partitions. Further, from (2.8) it follows that 


E {xi,x2, • • • = 53 • • • Dk, (2.10) 

9 

where 

= cum(xj,,xj,, • • • ,Xj„) p = 1,2,3, ••• ,g (2-11) 

where the ji are the elements in the partition Kj, and the summation in (2.10) is over 
all the possible partitions. Tables (2.1-2.3) show the moment-cumulant relations for 
orders 2. 3 and 4. The moment is the sum of all cumulant terms in the same column. 
Likewise the cumulant is the sum of all moment terms in its same column. The value 
of q and the corresponding partitions are also shown. 


partitions 

moment 

cumulant 

q 

Ki 

«2 


1 

1,2 



E{xiZ2} 

2 

1 

2 

Cli)(xi)Cl^)(*2)1 

-S{xi}£{x2} 


TABLE 2.1: Cumulants and moments for order k = 2. 


Let us denote a discrete random process defined on the integers — oo < n < oo 
as {a'(n)} or simply as x{n) when there is no possibility of confusion. Lnless otherwise 


11 



partitions 

moment 

Mif) 



Ki 

«2 

Kz 

D 

E!S9 



€^(* 1 * 2 * 3 ) 


D 


2,3 


C(^^(Xi)C^*1(452*3) 

-5{xi}5{x2X3} 

D 

2 

IB 


C(^^(x2)C^*H*1*3) 


B 

3 

IB 


CI^>(X3)C^"^(X1X2) “ 

-€{xz}E{xiX2} 

B 

1 


3 

C(l)(Xi)Cl^^(X2)Ct^^(X3)' 

£{xi}B{x2)£{xz} 


TABLE 2.2: Cumulants and moments for order ^: = 3. 


partitions 

moment 

cumulant 

q 

Ki 

K 2 

^3 

/C 4 

1 

1 , 2 , 3,4 




C(^l(XiX2X3X4) 


2 

1,2 

3,4 



QS^\xiX2)C^‘^\xzx^) 

-5{xiX2}5{x3X4} 

2 

1,3 

2,4 



C^^1(iiX3)C(^^(x2X4) 

-f{xjX3}5{x2X4} 

2 

1,4 

2,3 



C^^^(xiX4)C^^^(x2X3) 

— £’{XiX4}£^{x2X3} 

2 

1 

2,3,4 



C^^l(xi)C(^)(x2l3*4) 

-5{xi}5{x2X3l4} 

2 

2 

1,3,4 



C(^^(X2)C(®^(X1X3X4) 

-5{x2}£^{XiX3X4} 

2 

3 

1,2,4 



C(^l(X3)C^®^(Xil2X4) 

-£:{x3}5{xiX2X4} 

2 

4 

1,2,3 



C^^1(x4)C^®^(xiX2X3) 

-f{x4}f{xiX2X3} 

3 

1,2 ^ 

3 

4 


C(2)(X1I2)C^'^(*3)CI'^(X4) 

— £’{XiX2}^{x3}^{x4} 

3 

1,3 

2 

4 


C^*1(XiX3)C(^1(x2)C^^1(x4) 

2£{xiXz}S{xz)E{xa} 

3 

1,4 

2 

3 


C(2)(xa*4)C^i^(x2)C^^)(x3) 

25{XiX4}5{x2K{*4} 

3 

2,3 

1 

4 


C^^l(X2X3)C^^l(xi)C^^l(x3) 

2f{x2X3}f{xi}f{x4} 

3 

2,4 

1 

3 


C(2^(X2X4)C^'>(Xi)d*)(x3) 

2f{x2X4}£’{ll}5{x3} 

3 

3,4 

1 

2 


C(2)(X3X4)CI'^(X1)CI'^(X2) 

2£^{i3X4}5{xi}£:{x2} 

4 

1 

2 

3 

4 

C(1)(X1)C^'^(X2)C('^(X3)C^'^(X4) 

—65{xi}5{x2}£^{x3}f{l4} 


TABLE 2.3: Cumulants and moments for order A: = 4. 
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specified, the random process will be assumed to take on real values, i.e., x(n) € H. 
A random process {a:(n)} is defined to be m*^-order stationary if for each k < m the 
cumulant of order k for samples of the random process x(no),x(ni), • • • ,x(n*_i) is a 
function of only the lag differences = m — uo, h — ni — no,’", h-i = — no- 

For an m*^-order stationary random process, its A:*^-order moment and cumulant 
functions are denoted by 

= S{{Tin),x{n + /i),x(n + /2),-t- /fc_i))} 

( 2 . 12 ) 

k < m 

and 

Cx\h-l 2 ^-"Jk-i) = cum(T(n), r(n + + 4). • • ., 1(0 + 4_i)) 

(2.13) 

k < m 

For brevity, we will also use the notation and C^x\l) where / is the vector 

argument 

The size of the vector I is implicit with the order of the moment or the cumulant. 

When there is no ambiguity it is common to refer to the moment and cumu¬ 
lant functions as simply “moments" or “cumulants" of the process. For a zero-mean 
random process the second, third and fourth-order cumulants are given by 
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£ {x(n)i(n + /)} = (a) 

£ {i(n)x(n + li)i(n + /j)} = W 

£ {x(n)x{n + /i)x(n + /2)x(n + /a)} - - k) 

-cW(;2)c<’>(/2 - h) - c(“>((2)cf )(;■ - ' 2 ) 

(j. h) - - ( 3 ) - /<Si'>('3);3(f>(i2 - /,) 

-f-?>('3)/<S?’('. - ' 3 ) (c) 

(2.14) 

These formulas follow directly from the corresponding formulas in Tables 2.1-2.3 for 
random variables. 

2. Properties of the Cumulant 

Cumulants have properties that follow from the cumulant definition and 
are important for their application in digital signal processing and system theory. 
Among these properties are the following [24]; 

1. If x{n) is a Gaussian random process its cumulants of orders higher than second 
are zero. Since the moment generating function of jointly Gaussian random 
variables with mean vector and covariance matrix C is given by 

<^;r(u) = (2.15) 

the natural logarithm of this function (the cumulant generating function) can 
be considered as a power series with zero coefficients for terms corresponding 
to powers greater than 2. Thus the cumulants of all orders greater than second 
are zero. This is a property of any Gaussian process, white or colored. 


c «(0 = 
C«((„(2) = 
citKhJiM = 
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2. For a zero-mean random process, the third and fourth-order cumulants can be 
defined as 

= - £ {(j(n + li}g(n + It)- "gin + /j-i)} 


for /: = 3 or 4 


(2.16) 


where g(n) is a Gaussian process with the same second-order statistics as x(n) 
[23]. (This does not apply to cumulants of orders higher than fourth). There¬ 
fore, the cumulant can provide a measure of the process “difference from Gaus- 
sianity.” In general, if x(n) is a stationary random process with variance 
the coefficient of Kurtosis is given by: 

£{(x(n)-£{x(n)})*} 

( 7 ^ ^ 

_ cum(x(n),x(7?),x(n),x(n))-|- 3-cum(x(r7),x(n)) •cum(x(7r),x(7i)) 

__3 

_ c4%,0,0)-h3(Ci='^(0))2 


which equals zero for Gaussian processes because in this case C’x^^(0,0,0) = 0 
and Cx^^(O) = (7^. 

3. If a random process consisting of independent identically-distributed (i.i.d) ran¬ 
dom variables has a symmetric distribution (such as Laplace. Uniform. Gaus¬ 
sian and Bernoulli-Gaussian) its third-order cumulant is zero. This is not true 
for non-symmetrically distributed processes (such as exponential, Rayleigh and 
K-distributions). The coefficient of skewness is 


5{(x(n) £^{x(7?)}) I ^ cum(x(7?). x(n), x(r?)) C4^^(0,0) 


which is zero for symmetrically distributed processes. 


( 2 . 18 ) 



4. If each of the random variables x,-, i = 1,2, • • •, fc is scaled by a constant a,- the 
cumulant becomes 

cum(aiXi,a2X2,-"5<^fc^*) = cum(xi,X2,• • • (2.19) 

5. Cumulants are symmetric with respect to their arguments, i.e., 

cum(xi,X2,-*‘,Xik) = cum(xii,Xi,,--',Xu) (2.20) 

where (ii, ^ 2 ^ • * • i h) is any permutation of the indices {1,2, • • •, fc}. This implies 
for example that the third order cumulant of a stationary real random process 
has six regions of symmetry in the /i, h plane such that 

c(®)(/i,/2) = cithh.h) = c^^H-hJ2 - h) = 

- h--h) = - h) = - h.-h) (2.21) 

The principal region for which /i ^ ^2 ^ 0 is called the non-redundant region 
of support of If the value of the cumulant at any point (h^h) in 

this region is known, the values of the cumulant at corresponding points in the 
five other regions of support defined by (2.21) are also known (see Fig.2.1). 
The number of regions of symmetry increases rapidly with the order of the 
cumulant. The fourth order cumulant of a real random process has 24 regions 
of symmetry in the space of li.hJa - 

6. Cumulants are additive with respect to their arguments, i.e., 

cum(xo + Xj, X 2 . • • •, x*) = cum(xo, X 2 . • • •, x*) + cum(xi, X 2 , • • •, Xfc) (2.22) 

Therefore the cumulants of sums of random variables are equal to the sums of 
their cumulants. 


16 






7 . If constants are added to the values of the random variables the cumulant does 
not change its value, i.e., 

cum(ai + Xi, 02 + 3J25 •' • 5 Ofc + x*) = cum(xi, X 2 , • • •, x*) (2.23) 

where Oi, 02 , • • •,o* are constants. 

8. If the sets of random variables {xj} and {y^} are statistically independent then 

cum(xi + yi,X 2 + y 2 ,-'-,a;* + J/k) = cum(xi,X 2 , ••*,Xfe) + cum(yi,y 2 ,-"iyfe) 

(2.24) 

9. If a subset of the random variables {x^} for ?' = 1,2, • • •, A: is independent of the 
rest, then 

cum(xi,X2.---,Xfc) = 0 (2.25) 

This implies that for a sequence w{n) of independent identically-distributed 
(i.i.d.) random variables, the cumulant of any order is the multi-dimensional 
dtlta function, i.e., 

for w{n) i.i.d. (2.26) 

where is a real-valued constant and 6{l) is the unit sample function, i.e., 

f 1 for / = 0 
for//0 

Such a process is referred to as higher-order white noise process. 

3. Advantages of the Application of Higher-Order 
Cumulants 

Cumulants and their Fourier transforms (the polyspectra) have recently 
gained increasing interest in many diverse fields of application. However for decades. 
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due to the laek of analyticaJ and computational tools for the higher-order statistics 
applications, only the statistics up to the second order (correlation and power spec¬ 
trum) were used to any significant extent in signal and system analysis. Although 
the application of higher-order statistics involves an extensive amount of computation 
and requires a much larger set of data compared to the second-order approach, there 
exist applications for which the second-order methods cannot serve appropriately 
[23]. 

The second-order statistics are “blind” to the phase of a random process 
while the higher-order statistics reveal both amplitude and phase information. In the 
fields of linear system identification and process modeling, the solution of a problem 
is not unique if only second-order statistics are applied. The system that results 
from this approach can be either minimum-phase or one of many possible associated 
nonminimum-phase systems, because in any of these solutions the underlying second- 
order statistics are the same. When higher-order statistics are applied to the same 
problem the phase information that is revealed from the dat a can serve to determine 
which of these systems is the required solution. 

Cumulants, on the other hand, are “blind” to Gaussian processes. There¬ 
fore, cumulant applications are insensitive to any additive independent measurement 
noise if it is Gaussian. Thus cumulant-based methods can boost the signal-to-noise 
ratio when signals are (at least partially) corrupted with Gaussian noise. 

When the processes are non-Gaussian or nonlinear (i.e., generated from a 
nonlinear mechanism) higher-order statistics are essential in most cases. In real- 
world applications, processes which are truly non-Gaussian have traditionally been 
dealt with as though they were Gaussian due to the lack of appropriate mathematical 
tools. Similarly, nonlinear systems and models are frequently analyzed after being 
linearized in the vicinity of an operating point. In this case the solution is reliable only 
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within a range that ensures its stability (i.e., the region of convergence of the linear 
representation). The development of higher-order statistics tools, both analytical and 
computational, can and has provided solutions that are free from both Gaussianity 
and linearity limitations. 

When higher-order statistics are applied to problems in signal and system 
analysis, cumulants, instead of moments, are used for many mathematical and prac¬ 
tical reasons. Among these reasons are: 


1. For higher-order (non-Gaussian) white noise cumulants of any order (not mo¬ 
ments) are multi-dimensional delta functions (see (2.26)) and their Fourier 
transforms are multi-dimensional constant functions in the frequency domain. 
This has an enormous effect on the simplicity of the computation of system 
output statistics when the input is white noise. Moments do not share this 
property, in general. 


2. The cumulant of a sum of independent random processes is the sum of their 
cumulants (see (2.24)). The same is not true for moments. This enables working 
with cumulants as optraiors. 


Although the above applies to all higher-order cumulants, in many cases the 
third-order cumulant may not serve the purpose of the application and one may need 
to work with fourth-order cumulants. This is clear in the case of white symmetrically- 
distributed processes, for which the third-order cumulant is identically zero. The 
fourth-order cumulant is also used in cases of processes that have a relatively small 
third-order cumulant and a much larger fourth-order cumulant. 
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4. Moments of Jointly Gaussian Random Variables 


Let y = [yi, J/2, • • •, VhY be a vector of real-valued zero-mean jointly Gaus¬ 
sian random variables. Their joint probability density function is thus given by 

r (y) =-(2.28) 

(27r)t|Cl5 

where C is the covariance matrix with entries 

Cij = Cji = S {yiVj} = S {yjyi} (2-29) 

From the definition of the moment generating function (2.1) it is easy to show that 




(2.30) 


Expanding this in a power series yields 




m=0 


■lu^Cul" 




Tn=0 
oo 1 

= E 


k k 
jl=l j2=:l 


m=0 


m! 


1 \ m 
9 


oo 1 -2771 

= y —— 

777 ! 2 ^ 

m=0 


k k 

E E 

ii =172=1 

T m 

k k 

E 

ii=i 72=1 


(2.31) 


To find the value of the moment 5 {yiy 2 • • • y*}- the expansion (2.31) is compared to 
the Taylor series expansion of the moment generating function and the required value 
of the moment is found as the coefficient of the term j'‘uiU 2 • • • Ufc. This comparison 
reveals the following properties [24]: 

1. The sum of the exponents of the Uj's in (2.31) is always even and equals 2m, 
= 0.1.2, ••• . Therefore, in the Taylor series expansion, the coefficients 
are equal to zero if the sum of the exponents in u'^u'^ • • • Ufc* is odd or if the 
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number of variables k in u\U 2 "‘Uk is odd. Thus the moments of the product 
of jointly Gaussian random variables is zero if their number is odd. 


2 . If the number of random variables is even and equals 2p, where p is a positive 
integer, the required average of their product is the coefficient of U 1 U 2 • • • ujp in 
the term 

if** 1" 

(2.32) 


1 

pl2P 


> 1=1 


The required coefficient thus takes the form 


^ {vh vh'" y^p} = ^ Z n ^ {j/ji yh } 


(2.33) 


where the sum of products is such that such that none of the subscripts appears 
more than once in the product of the p pairs of the form £ {t/jiJ/jj} and the 
summation is done over all the possible pairing permutations. 


To simplify the expression (2.33) we first notice that £{yj^yjj} = £ This 

means that every 2’’ of the terms under the summation are equal and can be replaced 
by just one term multiplied by the factor 2^. Also, since the order of the terms 
in the product has no significance, there are p! identical terms for the same pairing 
configuration. As a result of these considerations, the required expectation of product 
of Gaussian variables in its simplest form is 


^ {yitth ■ ■ ■ yiri = 


(2.34) 


where the summation here is over all the possible distinct pairing permutations. The 
number of these distinct permutations equals 

(2p)! 


Aperm- 


(2.35) 


As discussed in Chapter III. the output of a “well-behaved"’ nonlinear system 
can be described by a series-like representation with or without memory. This series 
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can be considered as a multinomial of functionals of the input. When a moment 
(cumulant) of the output is computed, the average of the product of these functionals 
raised to some powers is encountered. As shown in Chapter IV, the output cumulant 
computation involves the expectation of the product of Gaussian variables raised to 
some powers in the form 


^ j (2.36) 


k 

The value of this moment is zero if the sum of the exponents i/j is odd. As 
explained in Appendix A, when the sum of the exponents in (A.l) is even the value 
of the moment can be computed by first constructing the k x k correlation matrix 


C = 


C(l,l) C(l,2) C(l,3) 
C(2,l) C(2,2) C(2,3) 
C(3,l) C(3,2) C(3,3) 


C(\,k) 

Ci2,k) 

C{3.k) 


(2.37) 


C{k,l) Cik.2) C{k.3) ••• C{k,k) 


where 


C{i.j) = SiyiVj} 


(2.38) 


and the same size matrix M of non-negativt integer entries such that 


M = 


2h 

?2 

is 

ik 


ffc+i 

2ik+2 

ik+3 

hk 


*2fc+l 

hk+2 

'2i2k+3 


(2.39) 

i(k-i)k+i 

i{k-l)k+2 

*(fc-l)fc+3 • • • 

2^2 ^ 



The matrix will be called the mvltiplicity matrix. Note that for clarity it is prefered 
to use k^ linearly indexed variables ij to represent the entries of M rather than 
variables with dual (column,row) index. We also force the diagonal entries to assume 
even values for reasons that are explained in Appendix A. The value of the moment 
is then given by 
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where the summation is over all the possible configurations of the matrix M such 
that the sum of its entries along the rows and the columns satisfy the condition 


i=l 3=1 

This means that the sum of the matrix entries in the row and column equals 
the exponent of yi in (A.l). The first product corresponds to the product of the 
diagonal entries of the correlation matrix C (the autocorrelation) while the second 
corresponds to that of the off-diagonal ones (the cross-correlation). 


B. POLYSPECTRA OF RANDOM PROCESS 

Moment and cumulant functions provide information about the different order 
correlations among the components of random processes. They demonstrate the de¬ 
pendence of these correlations on the values of the time lag differences among these 
components. In many applications information about the frequency content of the 
process, the power distribution over the frequency components, and the coupling 
between the different frequency components are of significant interest. This infor¬ 
mation can be obtained from the frequency functions known as the polyspectra. In 
this section we introduce the definitions of the polyspectra, their properties, and the 
spectral representation of a random process that leads to the procedure of computing 
the poly spectra. 
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1. Definitions and Properties of Polyspectra 


The multi-dimensional z-transform of the ^■‘^-order cumulant of a random 
process x(n) is defined as the Ar‘^-order complex cross-polyspectrum 

,Zk-\) = 53 ••• 53 ( 2 . 42 ) 

ll = —OO = —oo 

where and Zk-\ are complex-valued variables in the region of convergence 

of Sx^. The Fourier transform of the cumulant will be called the fc*^-order cross¬ 
polyspectrum and is obtained by letting Zi = namely ^ 

f; ••• E 

= —(X) —oo 

( 2 , 43 ) 

A sufficient condition for the existence of this quantity is 

£ ••• £ < oo (2.44) 

—OO —OO 

The polyspectra are defined for all real values of the frequency variables u-’j but are 
periodic in each variable with a period 2r. For the third-order case this quantity is 
called the bispectrum and for the fourth-order case it is called the trispectrum. As in 
the case of moments and cumulants we shall sometimes use the notation Sx ^{z) to 
denote S^\zi.Z 2 .’ • ■ ,Zk-\) and Sx\(^) to denote - where 


Z1.Z2. • • •, 


and 

U) = [a;i.u;2, • • • .a;fc_i]^ 

The properties of cumulants discussed in Section A.2. lead to corresponding 
properties of the polyspectra, namely: 

^The notation S^\ui,W2, -''rather than • •-.e^"'-’) b slightly abusive 

but will be used for notational simplicity. The meaning will be clear from the context. 
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Figure 2.2: Regions of symmetry of the Bispectrum of a real random process. 

1. The bispectrum, trispectrum and all higher-order polyspectra are zero for Gaus¬ 
sian processes. 

2. The bispectrum is zero for symmetrically-distributed processes which are i.i.d. 

3. When the random variables are scaled, their cross-polyspectrum is scaled by a 
factor that equals the product of the individual scaling factors. 

4. Polyspectra have symmetry properties similar to those of cumulants. For exam¬ 
ple, the bispectrum of a real random process has the twelve regions of symmetry 
shown in Fig. 2.2. 
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5. Polyspectra are additive with respect to their arguments 


c(fc) 

^(xo+Xl,*J, 


(2.45) 


6. Polyspectra are insensitive to the addition of constants to the random variables. 

7. The polyspectra of the sum of independent random variables equals the sum of 
their individual polyspectra. 

8. If a subset of a collection of the random variables is independent of the rest, 
their joint polyspectra are zero. 

9. For processes consisting of i.i.d. samples, the polyspectra are constant, i.e., 
independent of frequency. 

Similar properties also apply to the complex cross-polyspectrum S^\z) 

2. Spectral Representation of a Stationary Random 
Process 

The second-order statistics of a random process are fully described by the 
covariance function and its Fourier transform, the power spectral density function. 
These two quantities, are related by the Fourier transform pair 

L (2.46) 

1= —C» 

The power spectral distribution function Is a nondecreasing function of uj 

that equals zero at u.’ = -tt and is obtained by integrating the power spectral density 
function 




(2.47) 



The covariance function can then be described as the Stieltjes integral 

= ^ r (2-48) 

ZTT J^ir 

since 

= S^\u})du> (2.49) 

The representation is embodied in the “Wiener-Khintchine theorem [18]. The in¬ 
version of (2.48) gives 

=^=(0,+»)+f ^c(='{/)+f: (2.50) 

Similar to this representation there exists a spectral representation of the 
zero-mean random process x[n). The development of this representation starts with 
defining the complex-valued random process Zxi^) w'ith orthogonal increments such 
that 

e{dZx{^^)} = 0 (a) 

e{dZx{^\)dZxM} = (^) (2.51) 

where 

^fl foru.'i=u;2 (2.52) 

( 0 otherwise 

Then the zero-mean random process can be represented as 

x(n) = ;l r e^'^dZxi^) (2.53) 

ZTT J-T 

The interpretation of (2.53) is that the process x(n) is represented as a sum of com¬ 
plex exponentials in the continuous range [—tt.tt) with random amplitude \dZx{^)\ 
and random phase LdZxi^)- In this case, represents the mean-square am¬ 

plitude of the component of x[n) at frequency u;. The value of ^'x(^) is the average 
total process power and Sx[j^')d'^' is the contribution to the total process power from 
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the components in the frequency range between us and us + <iu;. The spectral repre¬ 
sentation of the random process is therefore related to the spectral representation of 
its covariance function (the power spectral density function). 

By inverting (2.43), the fc‘'‘-order cumulant is given by 


Jk-i) = 

7r)*“^ y~ir J-ir 

(2.54) 


(27r) 


As an example, the third-order cumulant for the zero-mean process is given by 
(2.14.b); so substituting (2.53) we have 


C^\k-,h) = £ {j(n)a:(n -t- li)x{n + k)} 

= -— r r r {dZx(u.'i)dZx(^'2)dZx(u.’3)} 

(27r)3 J-T J-T J-K 


(2.55) 


Since the process is assumed stationary, the cumulant value must be independent of 
n. Since this can only be true if u,’i -I- a ;2 + = 0; follows that 

C'i\kJ 2 ) = 7 ^. r r {dZx(a;2)dZ:,(^3)dZx(--'2 -u;3)} (2.56) 

\1TT y J — T 1 / —“TT 

Comparing this result with (2.54) we find that 


u,'2)dijJidu,'2 — S {dZx{'^i)dZx{us2)dZx{ *^ 2 )} 


(2.57) 


By a similar analysis it can be shown that the cross-bispectrum satisfies 


^(3) 

Xyw 


(u-'i . us2^duS\d‘jS2 — ^ {dZx('.<'’i )dZy(<x.’2)dZu,( u.'’i 


— u;2)} 


(2.58) 
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For the fourth-order cumulant we can write 

/2, U) = J;. h) - - U) 

- h) - - ' 2 ) 


— ^ gin(wi+c*^+W3+W4), 

(27r)^ /-ir •/-» /~ir i/~‘r 




X £: {(iZx (u;a )^fZx ( 0^2 MdZj, ( 0 ^ 4 )} 

- — r r 

(27r)* J-n J-ir 

(27r)2 y-ir J-ir 

-L_ r r {dZx(a;2)dZx(u.’4))} 

(27r)^ J-ir y-ir 

x-J— r r {dZx(a;i)dZx(u;3)} 

(27r)2 7-, y_T 

_L_ r r e"’‘(‘"*+"«)e^(">'’^5{dZx(u^3)dZx(u;4))} 

(27r)^ 7-ir y-ir 

x^ r r e^'”(“^+“^)e>“‘('>-'*)5{dZx(u.-i)dZx(u,’2)}(2.59) 
(27r)2 7_, y_,r 

In this case also to satisfj' the stationarity condition (2.59) becomes 




^ gi(wiJl+W2Jj+‘*'3is) 

(27r J—T «/—“TT j--‘K J—-K 

£ {dZx{^i)dZx{^2)dZx{^'3)dZx( *-*^2 ^’ 3 )} 


fJX 




£ {dZx(vtJi)dZx(—{dZx(a;2)dZx(—u.’2)} 6(u;2 -f u/'a) 


1 r r r , 




1 ^ 1 + ^*^2 ^ 2 + i*>3 ^ 3 ) 


e {dZ^MdZ^[-^2)} S {dZx(u;3)dZx(-u.'3)} <5(a;i -f cJa) 
^ f f f ^j{'^ih+>^2h+‘^ih) 

(2?:)® J-r J-n J-n 

£ {dZx(u.'3)(/Zx(—^’ 3 )} ^ {dZx(vt.'i)dZx(—i^’i)} -I- uJ2, 


(2.60) 


Comparing this result with the definition of the trispectrum by letting A' = 4 in (2.54) 
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yields the relation 


S^\uiTUJ2^^3)duJidu}2di*^3 = S {dZx{ijJi)dZx{ij‘}2)dZx{i*^3)dZx{—ijJi — u)2— 

—Z {dZx(wi)dZx(—wi)} Z {dZx(w2)dZx(—W 2 )} i>{y^2 + 

-£:{dZx(u;2MZx(-u;2)} Z {dZxMdZxi-i^z)} + t^a) 

-Z{dZxMdZx{-<^3)}S{dZxMdZx{-u;i)}6{u^+U2) (2.61) 

Similarly we can show that the cross-trispectrum of the processes x,y,w and v is 
given by 

Z {dZx{‘^i)dZy{u>2)dZfu{'jJ3)dZv{—u!i — uj2 ~ ^.’ 3 )} 

— Z'{dZx{i^'i)dZy{—u,'i)y6{uJi -j-iiJ2)Z {<iZu,(a’3)dZ„(—u-’s)} 

-Z{dZx(iOi)dZ^{-uJ^ )}<5(a.’i -f u^a{dZ„(-’2 )dZ, (-^’ 2 )} 

—£’{dZx(t<;i)dZ„(—a;i)} {dZy(u!2)dZy,{—Li!2)}6{u:2 +<^’ 3 ) (2.62) 

Now let us specialize these relations to the case of a zero-mean Gaussian 
process g(ji). Since for this process the cumuiants of orders higher than the second 
are identically zero, the left hand sides of (2..57) and (2.61) are identically zero. Thus 
for this process 

f {dZ„(^-i)dZ,(^'2)f/^p(-u.'i - u.’2)} = 0 (2.63) 

for all ui and u; 2 , and 

£ {dZp(u;i)dZg(u.'2)dZp(u.'3)dZg(—u-'i — o-’a — <*^ 3 )} = 

Z {dZp(u,'i)c?Zp(—o-’i)} ^(u;i -f- ^2)^{dZg{<jJz)dZg{—ij:2)) 

-\-Z {dZg{u!\)dZg{—U!i)}6{u!i U>3)Z {dZg{u.'2)dZg{—u!2)} 

+Z {dZp(a.'a )dZp( -^-1 )} Z {dZg{^2)dZg(-u>2)} (2.64) 



This means that for the Gaussian process the value of 

S {dZg{ux)dZg{uj2)dZgMdZ„{u}i)} 

is nonzero on the manifold defined by wi + 1*;2 + W 3 + W 4 = 0. Moreover on this 
manifold it is nonzero only on the three proper submanifolds defined by 

• (wi + W2 = 0) n (013 + 0)4 = 0) 

• (wi + u;3 = 0) n {u>2 + W4 = 0) 

• (u^i + u;4 = 0) n (012 += 0) 

The value is given by 

e {dZgiUl )dZg{0J2 )dZg{^>2)dZg{uJ^)} = 

5p(u;i)5ff(u;3)^(u,’i + u;2)6(a>3 + u.\)dijJiduJ2diCzdijj/i 
-|-Sp(uJi)5p(u;2)^(tfc’i + t^’3)^(^'2 + ^4t)dijJ\du}2d>jJzd'j:^ 

+5g (ii;i )S’j( 0 ^ 2 ) <5 ('^1 + <^'4)<^(<*’2 + u,’3)cL.’icL;2du;3du;4 (2.65) 

which has a three-dimensional region of support defined by u;i, uJz and u-’a; we use 
four frequency arguments for clarity however and the condition 

o'! -t- ^^2 + + <^’4 = 0 


to make the notation correct. 

Now assume we have a Gaussian process that is also white. Denote this 
process by io{n). In this case the right hand side of (2.65) becomes 

E {dZ,o{^\)dZ'ui{y^2)dZ.uj{i-^z)dZ^{>jJi)} = 

+ u;2)<5‘’(u.' 3 -f U’4) + + aJ3)^‘(i^’2 + t^’4) + < 5 ‘'(w^l + + t^s)] 

d^\ du,'2 du.'z da.'4 
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where is the continuous time unit impulse function, not to be confused with 

the dirac delta 8^,^ defined in (2.52) or the (discrete time) unit sample function 8{l) 
defined in (2.27). In general since the process u;(n) is zero-mean white Gaussian, all 
the cumulants except the second-order cumulant (the covariance) and the associated 
power spectral density function are equal to zero. Therefore when we apply this char¬ 
acteristic to any order polyspectra expressions similar to (2.57) and (2.61), equations 
(2.63) and (2.65) can be generalized as follows: 

{ 0 for odd k /« 

a^S^dii} for even k 

o Uf 

where we use the notation 

du) = dijJi djj\ • • • dLjJk 

and where the product is over the | pairs of and o-'j and the summation is over all 
possible pairing permutations. 

C. MOMENTS AND SPECTRA OF FUNCTIONALS OF 
A WHITE GAUSSIAN PROCESS 

1. Moments and Spectra of Linear Functionals 

For a zero-mean white Gaussian process U'(r7), the power spectral density 
function has the same value al for all frequencies; this is the average process 

power. Therefore the spectral representation of w is given by 

^ r (2-69) 

ZTT J—T 

where 

5{dZ„(u.-)} = 0 

e{dZ^MdZ*Jo:2)} = d'I'x(^'i)^«.u„ (2-70) 
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and 


= Sv,{uj)duj = aldu) (2-71) 

If this process drives a linear system with transfer function H{uj) then the output 
y{n) is another zero-mean Gaussian process that also has a spectral representation 
in the continuous region [-7r,7r) in the form of 

y[n) = ^ (2.72) 

The components of the input random process that comprise its spectral represen¬ 
tation are modified by the magnitude and phase of the transfer function to give 
the components of the output process with random magnitude |//(ju;)l|dZ„(a;)| and 
random phase lH{u;) -t- LdZ^(uj). In general 

dZy(u;) = H (2.73) 

Which satisfies conditions similar to (2.51) because 

S{dZy{^^)} = //(^■)^{dZ„(U.-)} = 0 

f{dz,(u.’i)dz;(a’2)} = /7(^i)//*(u’2)£:{dz„(.-i)dz:(^-2)} 

~ )l (2-^4) 

The output power spectral density function has the value 

5 (")(^’) = ( 2 . 75 ) 

If the system impulse response is denoted by h{n) then the output autocorrelation 
function is given by 

■Ry(0 = £{y{n)y{n + /)} = alrhil) (2.76) 

where r/,(/) is the system correlation sequence given by 

OO 

i^h{l)= hin + l)h{7i) (2.77) 
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which is the Fourier inverse of \H{oj)\^- 

If a set of linear systems with transfer functions ft(w),//r(u>),•••,/?»(“>») 
and corresponding impulse responses lii(n), Arin), • • •, ha(n) are driven y 
zero-mean white Gaussian process u,(n) with variance the outputs y.,ys.-.Va 
ate zero-mean Gaussian and have the spectral representations 

(2.78) 




; J 


= 1 , 2 , • • •, ^ 


for which 


dZy^ioj) = Hj{juj)dZU^) 

In this case the cross-correlation function is then given b> 


= ^{yii^)yM + 0 } - 


(2.79) 


(2.80) 


where 




n=-oo 


(2.81) 


The quantity ri,(/) will be called the syMem cross-cornlal.on sequence. The cross- 
Spectral density function is then gi\en bj 

= <^lsij{^')d'^ 


(2.82) 


where 


sij(^') = 4(^') = 

We refer to the last quantity as the system cross-spectral functwn. 


(2.83) 



2. Moments for Nonlinear Functionals of a White 
Gaussian Process 

A set of nonlinear processes Uj(n) can be formed by multiplying the outputs 
of the linear systems raised to some integer powers t/ij to give 

= yTi^)yTin) • • (2.84) 

The Uj are functionals of the white Gaussian process w{n) because each of the yj is 
given by 

2/i(”) = £ Wj{n-i)hj{i) (2.85) 

»= —OO 

The value of the m*^-order cross-moment function 

= €{ui{n)u2{n + li)--- Um{n -f Im-l)} 

= S{y?^ {n)y^Hn) ■ ■. (n + h)y^^n + I,)... y^'‘{nh) ■ - ■ 

4- lm-l)y 2 ”'H'n -f /m-l) + Im-l]} (2.86) 

is zero if the sum of the exponents is odd. For an even sum of exponents the procedure 
explained in the Subsection A.4. is followed; the mkx mk matrix is constructed 
as 



■R(0) 

R(/i) 

R(/2) 

•• mim-i) 


R^(/i) 

R(0) 

R(/2-/i) 

•• R(/„.-i-/i) 

CH) = 

u 

R’’(/2) 

R’’(/2-/i) 

R(0) 

•• R(/„_i-/2) (2.87) 


. R^(/n.-l) 

R^(/m-l - h) 

R^(/r„-l-/2) • 

•• R(0) 


where the k x k matrix R(n) is given by 


rii(n) ri2(n) ri3(n) ••• rifc(r7) 

r2i(n) r22{n) rjsfn) ••• r2fe(n) 

R(7?) = r3i(n) r32(r0 r33(n) ••• rzk{n) (2.88) 

rfci(u) rk2{n) rkzin) ••• rkk(n) 
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and rain) are the system cross-correlation sequences (see (2.81)). A multiplicity 
matrix M of size corresponding to is defined as 


M = 


M(/i) 

M{h) 

••• M(/„*-i) 

M(0) 

M(/2-/i) 

• • • M(/m-l ~ 

- h) 

M(0) 

• • • R(^m-1 “ 


M(0) 

M^(/i) 




M(0) 


(2.89) 


where each block M(/) represents the multiplicities of the r<,(n) in the block R(/). Al¬ 


though equation (A.15) is still applicable, we can benefit from the fact that C{iJ) = 
C{j, i) to save computations in this case and in similar cases where the matrix Cj, ^ 
is symmetric. The multiplicity matrix M can be modified because the two entries 
M{iJ) and M{j,i) result in pairing M{i,j) + M{j.i) elements from both yi and y,- 
to give Therefore we let both entries of M be equal such that 

each of them equals the value that their sum would have if the notation in (A.15) 
is used. The diagonal elements are not changed because they still have the same 
meaning. Subsequently the off-diagonal elements need to be examined in the upper 
triangle only, and the condition on the sum of the elements along the columns and 
the rows becomes that the sum of the tkmenis of the row equals the sum of the 

elements of the column and each equals i/j. 

Before proceeding, let us define some notation that will be useful in the 
remainder of this discussion. For a /:-dimensional vector v with elements 
and i/fc, we dfine the additive reduction or simply the reduction of v as the sum of its 
components. If we denote this by E[i/] we have the definition 

= (2.90) 

1=1 

For a kxl matrix M, with elements M{i,j) we define the reduction as a vector with 
elements equal to the sum of the rows of M. Thus 

i 

i;[M] = m where mi = ^M{i.j) (2.91) 
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Now since each of the cross-correlation functions equals the corresponding 
system cross-correlation sequence multiplied by the input variance <7*, and since the 
total number of the cross-correlation functions is always equal to half the total sum 
of exponents, then each term in the expression of the moment has the value (<t*) 
multiplied by the product of the system cross-correlation sequences rj that result from 
each pairing permutation. Therefore the moment in this case is given by substituting 
in (A.15) 


= S{ui{n)u2{n + /i) • • • Umin + Im-i)} 


»=i E[M]=i/ 2 » ij=i 

^ 11 hf{: ,■ \t 


mk 


mk 


{CtHjuj,}) 


Mjji til) 
7 




(2.92) 


a/0i.;2)! 

The condition S[M] = */ in the subscript of the summation means that the summa¬ 
tion is done over all the possible values of the matrix M that satisfy the condition. 


3. Spectra of Nonlinear Functionals of a White Gaussian 
Process 


In the frequency domain the (m — 1 )-dimensiona] moment cross-spectral 
function is defined as 


) 

oo oo 

= '£■■■ E tfrHh-h. 

with inverse transform relation 


m-1 


(2.9.3) 


>--/ 

m—1 




d^'id^'2 • • • d^'m-i (2.94) 
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Equating the right hand side of this equation and (2.86) yields 

TO —1 

dLj\du)2 • • • dwjji^x 

= f{ui(n)u2(n + /i) • • • Um{n + 

= W‘ WHn) ... yr‘(^)yr in + k)yrin + /i) • • • + /i) • • • 

+ + (2.95) 


Now substituting each y{n) with its spectral representation in (2.78) we obtain 
"" ' ’ 7 -" • • • dZyiu:,) 

U 

= i^T r "■ r • • ■ //(^v)e^"(“‘-'“'+-+“*') 

ZTT •/—IT •/—“JT 

V--- 

u> 

dZyj{u:i)dZ^{u!2) • • • dZ^{u;^,) (2.96) 

Let us define as new variables of integration that are necessary when we sub¬ 

stitute (2.96) into (2.95). Then the right hand side of (2.95) takes the form 


€< 


{:■ 


1 ^ir /’f / ”* * 

v> / •••/ n n n 


e.‘., Er;i7 ^i(E:., e:.-. e?.? 


where P equals the sum of all the exponents in (2.95). i.e. 


(2.97) 


P = E[»/] 


(2.98) 
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and must be even so that the expression is not equal to zero, 
condition implies that 

m k w,.. 


£ S Z) ‘^U**** = 0 


il=l t 2 =l 43=1 


The stationarity 

(2.99) 


Now interchanging the order of the integration and expectation in (2.97) produces 

(')- f ... r (n n 

ZTT J-v ./-ir / Uitatj J 

( 2 . 100 ) 


where the integration is taken over the manifold defined by (2.99) and the product 
notation used here is equivalent to the triple product expression in (2.97). By substi¬ 
tuting the expectation in the integrand with its value given by (2.67), the expression 
becomes 



/ t'l *-2 ^3 


( 2 . 101 ) 


Since 6^ makes this expression equal to zero outside the manifolds defined by the 
P!/((y)!2? ) permutations of the product of delta-dirac functions <S(u;i + a;, ); the ex¬ 
pression (2.101) is a sum of terms represented by f multiple integrations. Each term 
corresponds to a certain submanifold of that defined by (2.99) and it is multiplied by 
the number of permutations that result in the same integrand. The product of the P 
transfer functions H[^') is replaced by the product of y system spectral density func¬ 
tions 5ij(u.') and the sum of lag-frequency products /lU-v in the complex exponential 
is reduced correspondingly because 


r r 77i(M)7/2(ju;2)(5(^-i 

= r (2.102) 

Z 77 

When this result is applied to (2.101) we can develop a procedure to calculate the 
value of 5 ^"*^(^’i.u,’ 2. • • • ■ u-'m-i) given the linear system transfer functions (or the sys- 
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tern cross-spectral density functions) and the input variance 


4. Second-Order Cross-Spectral Function of Functionals 
of a White Gaussian Process 


Let us continue to consider the functionals Uj[it>(n)] of the zero-mean white 
Gaussian process u)(n) defined by (2.84) and (2.85). Note that in general these func¬ 
tionals do not have zero mean. The cross-correlation for any two of these functionals 

is given by 

/JHO = 1)} 

= + /)---yi^*(” + o} (2.103) 


and its frequency transform (the cross-power spectral density function) is defined as 


E (2.104) 

/ = —00 

If the sum of exponents in (2.103) has an even value P then the expectation is not 
identically zero and we can compute the value of the moment as described in (2.92). 
The right hand side of (2.95) is obtained as a sum of terms each corresponding to a 
possible permutation pairs of frequencies resulting in the expression (2.101). We start 
the procedure by constructing the 2k x 2k matrix of system cross-spectral functions 



S°(a;) S‘"(^’) ■ 

(S“)*^(^') S°(u.’) 


(2.105) 


where each term S°(t^’) or S‘^(u;) is of the form 

Sii(a,') Si2(^’) 513(0.') ••• 5ifc(a;) 

52i('^') 522(<^’) 523(u'') 52fc(o.') 

S(u;) = ■ 53 i('^’) 532(0;) 533(0.') ••• 53 fc(u;) ( 2 . 106 ) 

5 fcl(o.') 5 fe 2 (^') 5^3(0.') ••• 5 fcfc(o.') _ 
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with the system cross-spectral function Sij defined in (2.83). The significance of the 
superscript (0 or w) is explained below. A 2k x 2k multiplicity matrix M is also 
constructed in the same way as for the moment computation procedure. 

Let us now assume that the frequency pairing permutation in (2.101) re¬ 
sults in a component 5i^^‘(u;) of 5(^^(a;) which is obtained by summing all the possible 
components. In this component the frequency pairing is over all transfer functions 
that comprise the functionals uj and U 2 . Those transfer functions of ui do not have 
their frequency argument appearing in the lag-frequency product of the complex ex¬ 
ponential. Therefore if the frequency arguments of two such transfer functions are 
paired together they do not affect the lag-frequency exponent and do not have their 
system cross-spectral function multiplied by a complex exponential of its frequency 
argument. This spectral function is an entry of the upper left diagonal block S° of 
If the frequency arguments of two transfer functions that come from uj are 
paired together then the sum of these two arguments is zero and they therefore dis¬ 
appear from the lag-frequenc}' complex exponential. The resulting spectral function 
is an entry of the lower right diagonal block S°. On the other hand if a frequency 
argument of a transfer function of ui is paired with one of ii 2 , this frequency argu¬ 
ment appears in the lag-frequency product in the complex exponential. This spectral 
function is located in S‘^. 

For a specific configuration of the frequency pairing let us assume that pi 
system cross-spectral functions lie in the upper triangles of the two diagonal blocks 
denoted and p 2 are entries of such that pi -h pa ~ Then from (2.104) and 
( 2 . 101 ) 
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( —)a y ’ •' j ^Si{l>Ji)S2{uJ 2) • • • Spii^pi) 

P 

2 

xsp,+i(wp,+i)spj+2(wpi+2)-’-S|(‘^£)e'*^‘^'''^‘^ ^^dxJidu}2-’du}E 

(2.107) 

(Here we used a single digit in the subscript of the spectral functions for simplicity; 
they are actually indexed with two digits according to the indices of the correlated 
system outputs). Since each of the first pi spectral functions can be integrated 
separately, the result of this integration is a constant multiplying the rest of the 
integrals. The value of this constant is given by 


J_r 

27r J—K 


f Sij(u,')du: = rij{0) 

J —-K 


(2.108) 


w’here ru is the system cross-correlation sequence defined in (2.81). Therefore (2.107) 


reduces tc 


5(2)‘(u.>) = ((7^)2ri(0)r2(0)--TpJ0) e ' j__ 


+u ^pi +1; + 2^Pi 


■Spi+2(^>t+2)”--S£('^?)t 




(2.109) 


Now interchanging the order of the summation and integration we obtain 

= (<7o)"^i(0)^2(0)--Tp,(0)(^)»^ / ••• / 

Z/i tZ-'JT •/“‘TT 


^‘'p,-(-i(‘^pi-n)-''pi+2(-''i>i+2) • • • ('^’f) 2 ^ ^ ^ fl^Pi+1 • • ■ 

/= —oo 

( 2 . 110 ) 

where we have used the fact that formally the Fourier transform of results in 

an impulse 27r(5'(u,-) in frequency. The summation in the integrand equals 2» 6(a.’ — 
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- ajpi +2 • • • — wp) which reduces the number of integrations to p 2 - 1 over the 
manifold defined by 

Wpi+i+a;p,+2 H-|-a;£=u; (2.111) 

In other words we can conclude this development by the expression 

= (a^)?5®WconvW(5"(u;)) (2.112) 

where 5®^*^ represents the product of the system cross-correlation sequences (0) 
that are located in the diagonal blocks and included in the permutation. The 
value of conv(*)(5"(u;)) is obtained by performing multiple convolution on the system 
cross-spectral functions located in the off-diagonal block S" and included in the i 
permutation. Specifically, for p spectral functions Sj. S 2 . * • •, -Sp of S included in the 
permutation this operation is defined as 

convW(5‘"(c.’)) = 

(-^)p~^ f *•• / Si(u.’ — ^1 )-S2(^l ~ ^ 2 ) ■ ■ ■~ ^p-l ^p-l) 

2 ;: J-ir J-r 

V—---' 

p-1 

de^(I02 ■ - dOp-i ( 2 . 11 . 3 ) 


Finally, we obtain an expression similar to the expression of the value of the moment 


2 k 2* 

»=i SfM]=l/ 2 » ii=i 




(2.114) 


where u is the vector of exponents appearing in (2.103). 


5. Higher-Order Cross-Spectral Density Functions of 
Functionals of a White Gaussian Process 


By generalizing the procedure developed above, a procedure can be de¬ 
veloped to find an expression for the value of third, fourth and other higher-order 
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cross-spectral density functions. The third order cross-spectral density function is 
defined as 

= f; £ < 2 . 115 ) 

t 3=~00 

We start the procedure by constructing the 3k x 3fc matrix of system cross-spectral 


functions 


5?) = 


S® 

(S“>) 


•r 


S«i 

s® 


S«-j 

S«1- 

,T s° 


Uf2 


(2.116) 


where each of the k x k matrices is defined as in (2.106). The associated 3k x 3k 
multiplicity matrix M is also constructed. We assume that for a specific frequency 
argument permutation (specific configuration of the matrix M) in the expression 
(2.101) some system cross-spectral functions are entries of the diagonal blocks de¬ 
noted S®. The result of these functions is a constant value taken out of the inte¬ 
grations. The rest of the functions are divided into three groups. The first group 
consists of pi functions that are entries of the block . The second group consists 
of p 2 functions that are entries of the block S"*. while the third group consists of pa 
functions in the block Each of the associated frequency arguments appear 

in the lag-frequency product of the complex exponential. If the sum of the expo¬ 
nents of the Gaussian functions yi in the moment expression has an even value p, the 
corresponding component takes the form 
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\ t / fi = ~00 l2 = “<» 

£•••/: 511 (^^ 11 ) • • • -SlpiCt^lpi )52 i(W2i) • • •52p,(w2|.j)53l(‘^3l) * ‘ ’ ■S3p,(t^3ps) 

^ ~ —V-^ 

Pi -f P3 +P3 

.<i^3pj 

-L-L 5ll(wil) ■ • • -Slpi(u,’ipi )52 i(w2i) ' • • •S2pj(t^2p,)'S3l(‘^3l) ' * ’ '53pj(‘^3pj) 

>—-!✓ 

P1+P2+P3 

00 

^ (‘*'1-^11-Wipj “WJI- 

/l = -oo 

X ^ -"»P2+"«+•••+">« )cL>ii<fu;i2. duzp^ 

/2 = -CX> 

(2.117) 


where u;ij are the new variables of integration that are needed when we substitute 
(2.102) into (2.101). After substituting the two summations with the values (5'^(S[a;]), 
this expression equals zero outside the manifold defined by the two relations 


u^l — U.’ll 


“h • • • + + (cJ3i + • • • + u.’3p3 ) 


Uu\ — yjJ 


21 4--h a-’2p2 — + -1" ^'3p3 , 


(2.118) 


If we now define the argument 


0 = 1 x 131 + • • • + tUSpa 


(2.119) 


w’e have 


+ • • • + li’lp, — w*.’l — 0 


( 2 . 120 ) 


• • • + <x'2p, — x’2 + ^ 


( 2 . 121 ; 
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and the number of integrations is reduced to pi + P 2 + Ps ~ 2. This can be repeated 
successively as follows. The first group of functions is integrated pi — 1 times over 
the manifold defined by (2.120) to give 

27r J-ir J-T 

pl-1 

( 2 . 122 ) 

The second group is then integrated P 2 ~ 1 times over the manifold defined by (2.121) 
to give 

+«) = 

f ... / 52i(u; 2 + ^ — ^l)'S22(^l ■" ^2) • • ■'52p2(^P2--l)^^1^^2 • ’ • 

27r J-tr J-tr 

----' 

P2-1 

(2.123) 

and the third group is integrated on the manifold defined by (2.119) to give 

^1 )*^32(^l — ^2 ) * * ’ '^3p3(^P3-l )d^ld02 • • • dOp^^i 

(2.124) 



The last integral gives us the value of the required component 


=(<’=)? (n-s“<‘>) 

x(i) /' s®.- «).S® .,(-3 + 

ZTT J—ir 


(2.125) 


and the total value of the third-order cross-spectral density function is therefore 

3 fc 0 ( 3 ) 1 / 3 fc 1 3 fc 1 

5(V'..-3) = n*'(>)! S n 7«(3rir)T; H (2-i26i 

i=i r[M]=i/ 2"V^ ji=i ( 2 
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Finally an analogous procedure can be explained to compute the fourth- 
order cross-spectral density function 

5<‘>(u,.,u.,,W3)= E E E (2.127) 

/l = -00 I2 — — OO /s = -'00 

The Ak x Ak matrix of system cross-spectral functions is defined as 

'I 

ufi 

^•T 

and the same size multiplicity matrix is constructed to compute all the distinct 
permutations of frequency pairing and the number of their occurances. In this case 
we get six groups of the product of the system cross-spectra! functions located in 
the upper off-diagonal blocks. Let us assume that the number of spectral functions 
in the group is pj and their sum is P. The frequency arguments are denoted 
by In this case three intermediate frequency arguments ^ 1 ,^ 2,^3 

corresponding to the fourth, fifth and sixth groups located in g"*. S"'and 
Su>2-u>3 respectively are defined such that 

61 = u.’4i -f- tA,’42 -h • • • -h i*'’4p4 (2.129) 

O 2 — -f a;52 -f • ■ • -h u,’5p5 (2.130) 

63 = U>ei d" ‘^’62 "!-■••+ ^6pe (2.131) 

By procedures similar to those described for the third order polyspectrum the compo¬ 
nent 5 '^^^‘(a’i,a; 2 ,i^’ 3 ) is obtained by performing the following multiple convolutions: 


s(^) = 


S" S“ 

(S<-I).T go 

(S"=>)’^ (S 

(guijj.T (g 


gWi—W2 gWl— 

gO gt*;2‘“‘*^3 

^*7’ gO 


(2.128) 


48 



— ^2) = 

' j T -^1-^2- 1?l)5l2(t?l - t?2) • • • Sip, (t?p,_l)<ft?l<il?2 • • • <^t?p,-l 

p,-l 

^^v2{^2 + ^1 — Oz) = 

( —)*^ "'J S2i(u;2 + ^1 — ^3 —’^1)522(1^1 — ^^2) •'•S 2p,(l?p,-l)<^l?lrfl?2 ••‘t^l^pj-l 

'-V-^ 

Pi -1 

+ 02 + 03 ) = 

( —)^* •••J S3i{u;3 + O 2 + O 3 — ■0i)s32{'0i — Oz) •' • S3j^{'dp^^i)dOid-d2-••d'Op^^i 

P3-1 



f S4i(01 — t?l)S42(l?i — i92) • • 

— TT 

■ S4p4 (l?p4 —1 )<fl?i c/??2 ■ 

••d??P4_i 

P4-1 

sl^M) = 

1 

f S5i(02 — t^l)S52(l^l — 1 ^ 2 ) • • 

— IT 

• S5pb(^Ps- 1 )dt?ic/l?2 • 

••di?p^_i 

Pb-1 

^S..6(»S) = 

1 

^ S6i(03 — t?l)562(l^l — 1 ^ 2 ) • • 

— IT 

•■S6p6(^^P6-l)«^^lO'l?2- 

^^P6-l 


P 6-1 


(2.132) 


Finally the spectral function component 5^^)‘(a;i,u.’2,u;3) is obtained by performing 
the triple integration 


SWV„U...W3) = (<Tj)f(n5'>W) (;!' 


xS'L;(‘-'! + «. - «3)5S..3(->3 + «, + «3)52„<(«l)5‘2.»5(«2)-Si2- 






(2.1.33) 
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and the expression for the fourth-order cross-spectral function becomes 


Ak 


Si*\u;i,u;2,u;3) = JJ *^(0^ 

r[M]=i/ 


Si^u;) I* 


t=i 






4 k 2 


(2.134) 


In the same way cross-spectral density functions of any order for the functionals of a 
white Gaussian process can be obtained through a procedure of successive convolu¬ 
tions. 
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III. THE WIENER MODEL FOR A 
DISCRETE NONLINEAR PROCESS 


A discrete time random process can be modeled as the output of a discrete 
system driven by another discrete time random process whose statistics are known 
up to a sufficient order. This system is generally specified by a set of parameters that 
can be of either finite or infinite size. The model is successfully constructed when the 
system parameters are chosen such that the statistics of the system output match the 
statistics of the modeled process up to a required order. Here the underlying system, 
although nonlinear, is assumed to be stable and causal (non-anticipative). This 
assumption satisfies the realizability requirement of the model and (as shown later) 
generates a structured s 3 -stem architecture. In this chapter we present the Volterra 
series representation of discrete nonlinear systems and then develop the Wiener model 
of discrete nonlinear systems and their associated random processes. Although most 
of the published work on the Wiener model is formulated for continuous time, we 
adapt the theory in this dissertation to the discrete time case. Consequently, there 
are several extensions and new results that are developed in this chapter prior to the 
main results of this dissertation which are presented in the later chapters. 

A. THE VOLTERRA SERIES REPRESENTATION OF 
NONLINEAR SYSTEMS 

The output {x(7?)} of certain “well-behaved’’ nonlinear systems can be related 
to the input sequence {?c(r?)} by a Taylor series-like representation developed by the 
mathematician Vito Volterra and first applied to the study of nonlinear systems by 
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Norbert Wiener [22]. The details of this series representation of nonlinear sj'stems 
can be found in several places [18, 22, 17]. We give a brief outline of the theory here. 

1. Higher-Order Volterra Kernels and Operators 

The Volterra series representation for a causal nonlinear system with input 
ii)(n) and output x(n) can be written as 

oo oo oo 

x(n) = ho + hi{k)w{n — k) + ^2 ^2 k 2 )w{n — ki)w{n — ^ 2 ) 

*=0 *1=0 *2=0 

0000 

+-1- • •• hp{ki,. ..,kp)w{n - ki)---w{n - hp) + • • • (3.1) 

*1=0 *,=0 

where hp{ki^.. .,kp) is called the p‘^-order Volterra kernel which is identically zero 
for any of the ki < 0. An equivalent way of representing the system is by a series of 
operators 

x(??.) = ?fo[«’(”)] + Wi[i^’(”)] + '^ 2 [»'(”)] H-h 'Hp[w{n)] H- (3.2) 

where 

Hpluin)] XI ■ ■ ‘ .^>)u’(n - /ri) • • • u’)?? - kp) (3.3) 

*1=0 *,=o 

is the p*^-order Volterra operator or functional. In this respect the following is worth 
noting: 

• The Volterra series is a power series with memory. If the input is scaled by a 
factor a the new output is given by 

.t(7 ?) = ■Wo[u’(n)] + a?fi[u’(n)] + a^'H^luin)] +-h a^Hp[w{n)] H- ( 3 . 4 ) 

Thus it is inherently nonlinear. 

• There is no loss of generality if the Volterra kernels are considered to be sym¬ 
metric with respect to their arguments because, if they are not. there e.xists a 
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procedure by which any asymmetric kernel can be made symmetric. In par¬ 
ticular, for any asymmetric kernel .. ,kp) the symmetric one is given 

by 

.fc,) (3.5) 

where the summation is taken over all the p\ possible arrangements of the ki s . 
Any system constructed using asymmetric kernels can also be constructed with 
the symmetric one [17]. 


The individual Volterra functionals are homogeneous since 


'Hp[aw{n)] = a^7ip[w{n)] ( 3 . 6 ) 

Thai is, a change in the input magnitude results in a change of the output 
magnitude but not a change of the output waveform. 

To discuss the Volterra system in the frequency domain let A"(u;) and 
be the Fourier transforms of x{n) and u’{7j) respectively.^ Thus X(uj) and x{n) are 
related by the Fourier transform pair 

n=:~00 

— TT i X[u;)e^'^du) 

ZTT J —“jr 

and lT(u;) and w{n) are related similarly. Further, let the generalized transfer func¬ 
tions be defined as the multidimensional Fourier transforms of the kernels 

=£•••£ .( 3 . 7 ) 

fci=0 fc,=o 

^For the present discussion w(n) and a(n) are assumed to be deterministic and the corresponding 
Fourier transforms are assumed to exist. 
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which are symmetric functions of the frequency arguments Wi if the kernels hj, 
symmetric. The response of the order Volterra operator 


are 


Xp{n) = 'Hp[w(n)] 

be obtained by first defining the associated (artificial) functional 


oo oo 


* 1=0 k ,=0 


and its p-dimensional Fourier transform given by 
A^(p)(u)i,a>2, • • ’ ,‘^p) = Hp{u;i,u)2', • • • 


(3.10) 


Then the response of the p‘'‘-order Volterra operator is given in the frequency domain 
by [17, 18] 

A^p(a;) = _ f •••/ A'(p)(w — ^1. ^1 — ^2) • * * 1 (3-11) 

(zTTjP ^ J-ir J-'K 

which is a (p — l)-dimensional convolution, generally not simple to evaluate. The 
system output in the frequency domain is then given as the sum of the terms 


A^(.^) = Ao(-') + A'i(a,') + • • • + Ap(w) + • • • 

2. Limitations of the Volterra Representation of 
Nonlinear Systems 


(3.12) 


The practical application of Volterra series in nonlinear system theory has 
two main difficulties. One concerns the measurements of the Volterra kernels and 
the second concerns the convergence of the Volterra series. 

The measurement of a general system’s Volterra kernel is possible only if 
the contribution of each of the \olterra operators can be separated from the total 
system response. For the work done in this direction, the nonlinear system to be 
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identified is assumed to be representable by a finite order Volterra series and each 
kernel is obtained by isolating the higher order operators. There is no exact method 
to isolate an individual Volterra operator for a system that is not of finite order unless 
approximations are assumed possible. 

The problem of convergence of the Volterra series representation of a non¬ 
linear system is similar to that encountered in the Taylor series representation of a 
function. This is evident from equation (3.4). This representation is expected to 
converge only for a limited range of the magnitude of the system input. Moreover, 
nonlinear memoryless systems that include saturating elements cannot be character¬ 
ized by a Volterra series that converges for all magnitudes of the input [17]. 

The limitations of the power series representation of functions are circum¬ 
vented by using orthogonal functions. In power series representation of functions the 
value of the series expansion of the function is strictly convergent for the values of 
the argument; therefore a region of convergence must be defined. In orthogonal func¬ 
tion representation, the functions are required to converge in the mean square sense; 
therefore the convergence requirement is not imposed on the argument value but on 
the average of a product of functions. Wiener extended the notion of orthogonality 
to functionals and built a canonical model for nonlinear systems. This extension, 
known as the Wiener model, is discussed next. 


B. THE WIENER MODEL FOR A DISCRETE TIME 
NONLINEAR SYSTEM 

Weiner defined a class of systems suitable for his representation w’hich has been 
called the Wiener class. A system that is a member of the Wiener class has the 
following properties: 
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1. The system is not “explosive.” For a finite variance random input the system 
output also has finite variance. 

2. The system does not have infinite memory. The present value of the system 
output must become asymptotically independent of the past values of the input. 

1. 6-Functional Representation of Nonlinear Systems 


From the Volterra functionals, Wiener formed a set of orthogonal function¬ 
als which he called the G-functionals, since they are orthogonal if the input is a 
white Gaussian time function [22, 17]. The type of convergence of the orthogonal 
series is convergence in the mean. As a result, the class of nonlinear systems that 
can be represented by the Wiener G-functionals is larger than the class of systems 
describable by a Volterra series. Although Wiener developed and reported his results 
for continuous time, we shall state the results in discrete time as is more suitable to 
our purpose. Let us define the functional 


rp[/?p. /!.p_i(p), /;p_2(p). • • •. ho(p). tr(77)] 

p 

^p) + 


* 1=0 fci =0 


ki)w{n — ki)w{n — /:2) • • • w{n — A\)J 

(3.13) 


as the degree non-homogeneous \'olterra functional. Observe that Fp is a sum of 
homogeneous Volterra functionals of the type defined by (3.3). The p in parenthesis in 
the subscript of a kernel indicates that it is a member of a degree non-homogeneous 
functional, and we define /ip(p) = hp. Wiener developed the set of the orthogonal G- 
functionals by requiring that^ 

^ {7i,nN(n)]rp[/ip./!p_i(p)-V 2 (p)----'^o(p):«i’(n)]} = 0 for m<p (3.14) 

^Although the original results were formulated using temporal averages we will restate the results 
in terms of ensemble averages. 
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when the input w{n) is white Gaussian noise with variance al- In other words, Fp is 
orthogonal to any functional Tim of lower degree. The particular functionals Fp sat¬ 
isfying these conditions are called G-functionals and are denoted by Gpigp, u'(n)]. The 
quantity gp is the leading kernel of the functional from which the kernels 9p-i{p),9p-2(p)-, 

* ’ • i5o(p) are derived. The symmetry properties of both the Volterra and the Wiener 
kernels and the average of the product of Gaussian processes result in the following 
properties: 

1. For the even order G-functionals, the derived kernels of odd order are zero and 
the functional expression in the form of (3.13) has only even order kernels. 
An analogous statement applies to the odd order G-functionals for which the 
even order derived kernels are zero. The general expression for the Wiener 
G-functional is then given by 

[EJ 

Gplgp-.wiii)] = Gp_2mN(n)] (3.15) 

m=0 

where [|J represents the greatest integer less than or equal to | and Gp_2m[uH»^)] 
is the (p — 2mY^ degree homogeneous functional given by 

Gp-2m[u’(n)] = 
oo oo 

Y S 9p-2m{p)(h- h- • • •. V2Tn)u’(n - ki)w{n - kz) • • • w{n - kp.2m) 

(3.16) 


2. The kernels Pp_ 2 m(p) for m = 0.1,2,..., [f J are derived from the leading kernel 
gp by the relation 




{-irp\{air 


(p-2m)!(m)!2”‘ 


^ ^ * * ‘ ^ ^ gpi^l • • ' ' ' •! •. k\, k2, ‘ ‘ ‘ , kp-2m) 



(.3.17) 



The coefficients of the integrals of the derived kernels of Cp are similar to the 
coefficients in the expression of the Hermite polynomial 


H,(i) 


If (-irp!(tr;r /p-r-i 

^ (p - 2m)!(m)!2™ 


(3.18) 


Both of these properties are shown in Appendix B. 


In general, the Volterra functionals Hp and the Wiener functionals Qp of the 
same system are different. The Volterra representation of a general nonlinear system 


is 

l(„) = f;«iMn)] ( 3 . 19 ) 

»=:0 

while the Wiener representation of the same system is 

(3.20) 

i=0 

(In both cases the sum is possibly finite.) Since the Volterra functionals are homo¬ 
geneous, the Volterra functional of some degree p in (3.19) is obtained by combining 
all the p*'* degree kernels in (3.20). For a finite order system the Volterra functional 
is obtained by summing the leading and the derived functionals of the same degree 
along the columns of Table 1.1 while the G-functionals are formed by summing along 
the rows in the same table. The G-functional expansion of "Hp is therefore 

OO 

'Hpluin)] = Y. Gp{p+2m)[9p{p+2m)\v:(n)] (3.21) 

m=0 

and the kernels hp and Pp(p+ 2 m) are related in an analogous way. 

From the above considerations it is important to observe the following: 

1. The Wiener functionals and kernels are functions of the power level of the input 
(tI w’hile the Volterra functionals and kernels are invariant with the power level 
of the system input. Further, the Wiener G-functionals are Mi orthogonal 
unless the input is white Gaussian with variance 
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G functional 


homogeneous functionals 

Go 

Go 


Qi 


Gi 

Q 2 

Go(2) 

G 2 

G 2 


Gi(3) G 3 

Ga 

Go(4) 

G2(4) ^4 

Gi 


Gi(5) G3(5) 


Volterra fuctionals: "Hq "Hi 'H .2 ^3 ^4 ^5 

TABLE 3.1: The relation between the Wiener and Volterra functionals 

2. The class of systems describable by the Volterra representation is a proper 
subset of the class of systems describable by the Wiener G-functionals. Since 
the Wiener G-functional series is an orthogonal series for the white Gaussian 
input, it does not require the restriction on the magnitude of the input to 
avoid the strict convergence problem as in the case of Volterra series. This 
means that a system may have a G-functional representation and not a Volterra 
series representation. Any system that can be represented by a Volterra series, 
however, can also have a G-functional representation. 

Before continuing the development of the Wiener theory, it is necessary to 
point out two important additional properties of the G-functionals. First, it follows 
from the general form of the functional and (3.14)-(3.17) that the G-functional is 
linear with respect to its kerneL i.e.. 

Qp\cigp + C 2 /p: u’(n)] = c-iQp\gj,\ w{n)] -1- czGplfp. u.’(n)] (3.22) 

where Ci and C 2 are any real or complex constants. Secondly, the expectation of the 
product of any two G-functionals (say for two different nonlinear systems) is given 

‘See [17] for further discussion. 
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by 

S{gpA9pr'^Hn)]Qp,[fp,-.Mn)]} = I for i ^ (3.23) 

where 

7pi(*^o) = " ’ 5^ 9pi{ki,k2^’- ‘ ikp^)fpi{ki,k2.,'• • ikpi) (3-24) 

* 1=0 <^ 1=0 

The latter is equivalent to taking the expectation of the product of the two leading 
kernels only. The proof of these statements is also given in Appendix B. 

2. Orthogonal Development of the G-Functional 


It is shown (in Section C of this chapter) that the \\ iener kernel Qp can be 
expanded in the p-dimensional discrete Laguerre series 


oo OO 

^ 2 ^ • • • ? ^p) — *’* S (^l)^m 2 (^ 2 ) • ‘ ■ ^mp(^p) 

1711=0 mp=0 


(3.25) 


where \m{k) is the m*'* degree discrete Laguerre Junction and that the G-functional 
of can be correspondingly expanded in the form 


for white 
by sums 


LfJ 




(-i)V(^, 


^ Xi(k)w{n - k] 


(p-2m) 


(3.26) 


ap[V;u>(n)] 

noise with variance cr*. It is thus seen that gp[gp' i^’(^)] can be represented 
of products of terms 


yi(n) = 'Z^Ak)uin-k) (3.27) 

fc =0 

where yi{n) is the output of a system whose impulse response is the degree La¬ 
guerre function Xi{n) and whose input is u'(ri). This output is also a zero-mean 
Gaussian process with variance <7^. Thus from (3.26) it is seen that Qp can be rep¬ 
resented by a set of linear systems that provide all of the “memory” followed by a 
nonlinear system u'ithout memory that acts to form sums of products of the yiin). 
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section 1 section 2 section 3 



Figure 3.1: The Wiener model of nonlinear processes. 

Further, the coefficients in (3.26) can be identified as those of the Hermite polynomi¬ 
als. Therefore (using the linearity property (3.22)) the G-functional is simply given 
by the equation 

Gp[ 9 p\uin)] = 2 ••• Cmi-m, n Hp,(ym^(r?)) (3.28) 

mi =0 mp=0 j=l 

in which m j^ ^ mjj for 71 / 72 and Pj — P- The Hermite polynomials are 
orthogonal and the expectation of the product of two polynomials is given by 

£{Hp,(j/„,(n))Hp,(?/„,(n))} = (3.29) 

Since the Wiener representation of a nonlinear system is the sum of G-functionals, 
it follows that the system has the same form of representation as the individual 
tr(??)]. If the nonlinear system model is driven by the required white Gaus¬ 
sian noise, we have a rich representation for a large class of random processes. Fig. 
3.1 thus represents what will be referred to as Wiener model of nonlinear processes. 
This model consists of three sections. Section 1 is a single-input multi-output linear 






system with memory that can be represented as a bank of linear filters driven by 
the same input w{n). The input w{n) is a zero-mean white Gaussian process with 
variance cr^. The impulse responses of the filters in section 1 are the discrete Laguerre 
functions; 

Fig. 3.2 shows the structure of this section which, in general, has infinite length. 
Section 2 is a multi-input multi-output system with no memory, which can be repre- 



Figure 3.2: Section 1 of the Wiener model 

sented by an infinite bank of identical memoryless nonlinear single-input multi-output 
blocks, each driven by one of the outputs of section 1. Each block of this section con¬ 
sists of parallel nonlinear functions described by the Hermite polynomials. Although 







in the general model this section may also be infinite, processes with a finite degree 
of nonlinearity will require sections only up to some finite order Nff- The structure 
of one of the nonlinear blocks is shown in Fig. 3.3. Section 3 is a multi-input 



Figure 3.3: Section 2 of the general Wiener model. 

single-output memoryless system whose inputs are the outputs of section 2. The 
model output is the weighted sum of certain products of the inputs. The products 
are formed by taking one and only one output from each block in section 2. Sec¬ 
tions 1 and 2 are the same for all systems and therefore all processes in the Wiener 
class. The parameters, which are the weighting coefficients of section 3, characterize 
a particular process. In this case a stationary finite-moment random process can be 
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approximated by a finite size Wiener model [17]. The model size is given by : 

1. The highest order, Nl, of the Laguerre functions in section 1. This section con¬ 
sists of the iVt-f 1 systems with impulse responses given by Ao(n), Ai(n), • • •, 

It will be seen shortly that the order Nl roughly corresponds to the “length of 
incinory” of the nonlines-r system. 

2. The highest order, Nn, of the Hermite polynomials in each of the nonlinear 
blocks of section 2. Therefore, each block consists of iViv -h 1 parallel nonlinear 
memoryless systems described by the Hermite polynomials Ho(yi(f)), Hi(yi(<)), 

• • • As mentioned earlier, this corresponds to the degree of non¬ 

linearity present in the process (quadratic, cubic, etc. ). 

Let us now consider a particular vector of indices of size Nl + 1 defined as 

ai = [Qio.Oii,---,aiArJ (3-30) 

with 

aij € {0,1. • • •, Nn} for j = 0.1, • • •, A'l (3-31) 

Then the model output is given by 

x(n)= E (3.32) 

i[ai]=o 3=0 

in which the summation is over all vectors a. such that the sum of the components 
is less than or equal to A'jv- For a particular ar the term 

j =:0 

is a multidimensional Hermite polynomial in the outputs of the linear section. We 
can refer to this multidimensional Hermite polynomial as the Q-polynomial Thus 
the model output can be written more concisely as 

x(n)= E ca,Qa,(n) (3-34) 

i:[ai]=o 
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The quantity S[ai] is the order of the functional represented by the Q-polynomial 
Qori(^)- The Q-polynomials are zero mean except for the one with index ao = 0 
where the mean equals one. Then 

S {x(n)} = Co (3.35) 

C. DEVELOPMENT OF THE LINEAR SECTION OF THE 
DISCRETE WIENER MODEL 

A discrete nonlinear system that is a member of the Wiener class of nonlinear 
sytems can be represented using the G-functional representation. The kernels of the 
G-functional have the property 

E (3.36) 

ki =0 kj —0 kp=0 

Therefore these kernels can be expanded in terms of a complete set of orthonor¬ 
mal functions. In this section we derive the expression for these functions which 
were shown to be the transfer functions of the bank of the linear systems that rep¬ 
resent section 1 of the model. We also develop the orthogonal expansion of the 
G-functionals. 


1, The Transfer Functions of the Linear Bank, the 
Discrete Laguerre Functions 

The single-sided discrete time function Vm(?^) of order m is defined as 


Vm(«)= 


7 ?! 


-p" "*u(n - m; 


IH<1 


(3.37; 


W 


rrt!(ti — m)l 

here u(? 7 ) is the unit step function. The ^-transform of easily shown to be 


(z - 


for |r| > \p\ 


(3.38) 
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Since the functions v^in) can be thought of as individual terms of a polynomial, it is 
resonable to expect that an arbitrary time function, say <t>{n), could be represented 
in some way by a weighted sum of the The problem is more tractable however 

if we instead use a set of orthonormal functions, derived from the Vjn{'n), in the 
representation. These functions denoted by Aj(n) can be derived as follows. 

The functions Ai(n) are formed by a linear combination 


At(n) — Cn»tt'Tn(^) 


(3.39) 


with 2 -transform 


Aj(2) — CnitTm(~) 


(3.40) 


The coefficients Cmi for m = 0,1, • • •, t and i — 0,1,2,3, • • • need to be determined 


such that 


which is equivalent to 


^ Ai,(r7)Aij(77) = Si 


i z K\i,iz).\i,{z '^)dz = Si,i, 

ZTTJ JCI 

where the contour Cl is the unit circle in the complex 2 -plane. 


(3.41) 


(3.42) 


To satisfy the conditions, we proceed as follows. Substituting (3.38) in 


(3.40) we obtain 


_ Coi2(2 - pY -f Cij2(2 — p)*~^ -f-1- Cj,; 

(--/>)‘+^ 

zPd = ) 


= Coi 


- />)*' 


where Pi{z) is an degree polynomial in 2 that can be factored as 

P,{z) = {z-Z^){z-Z2)---{Z-Zi) 


(3.43) 
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Then 


K{z) = Coi 


z(z - Zl)(z - 22 ) ■ ’ ' (z - Zi) 


(3.45) 


--.V-/ -in 

is a rational polynomial function which has one zero at the origin, i zeros at 21 , 22 , • • •, 2 j, 
and (i + 1) poles at p. Now since 






(3.46) 


it follows that 


27rj Jci 


‘Ai,(2)Ai,(2 ^)C?2 = -Coi,Coi,/0' 


J_ f (2 - 2i)(2 - 22) • • • (2 - 2,,)(2 - 2r^)(2 - 22-^) • • • (2 - 2^^ 


i7r_; Jci 


(2 — p)’i+^(2 — 


(3.47) 


which in general is not zero. If the coefficients c^a. cu-• ■ ■. ca in (3.43) are chosen 
however such that the zeros 


Z\ — Z 2 = ••• = p 


then (3.45) becomes 


Ai(2) = Coi- 


■[z-p 

{z-pf 


I--I > 


(3.48) 


and (3.47) becomes 


27rj\i/ 'Ai.(2)Ai,(2-')d2- 


, 1 / (2-/)-^)»>-»»- ^ 
2rji fci {z — 


"QHiCOtjP 


-dz (3.49) 


If ?'i < ?'2 the integrand does not have any poles inside the unit circle, which means 
that the integral equals to zero. If zi > ?2 then it has 2 i -22 + 1 poles at p but the 
integral is also equal to zero since 


1 

[21 - 22)! c? 2 ‘ 4-‘4 


.-1 \tl-Xj-l 
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Therefore the functions Ai(n) whose transform is given by (3.48) satisfy the orthog¬ 
onality condition. 

To normalize this set of orthogonal functions, set ii = h = i m (3-49) and 


equate to 1: 

z-'\i{z)Uz-')d 

2x7 Jci 


Z = -CoiP 


I -- -dz 


2xj Jci {z - p ^){z - p) 




Thus it follows that 


= 1 



(3.50) 

Coi = ^’V 

— p^ 


(3.51) 

P P (z - 

■p-^Y 

p)‘+i 

kl > 1^1 

(3..52) 


= VI - P^P ^ 

To obtain an expression for Cmi in (3.39). equation (3.52) is put in the form 


Ai(-) = Vi-^v 


{z - pY 


= 


= v'l-z’V 


—y(_i)>-'"(! -p^ 

_ rtV+i ^ 


* . 7m+l 

m=0 


- V ^ _ p).+i • V - 

* «, 2mil f i \ ^ 

= ^ ("I) P (1“^) ^ p)m+l 

tnrrO 

= E(-i)V-"‘(i-^^)^(L 

m=0 ^ ' 

Equating the coefficients of Vml^) in (3.40) and (3.53) we obtain 


{z-pY 


(3.53) 


= (-1)V’"(1 -p 


3m±J / ? 


(3.54) 


m 



Then substituting (3.37) and (3.54) in (3.39) we arrive at the explicit expression 
Ai(n) = - "■> 

= \A^ E(-ir (t,) (^) (1 - - ™) (3-55) 

This is the discrete Laguerre function \i{n) which describes the impulse response of 
the system in the bank of linear filters in the discrete Wiener model. The bank 
of linear filters comprising section 1 of the model is then represented by the impulse 
responses described by the discrete Laguerre functions. 

Basic results on the discrete Laguerre functions can be found in [25, 26, 27]. 
This complete set of orthonormal functions has recently gained considerable interest 
in linear system applications. Among these applications are the problem of iden¬ 
tification of causal stable systems [28. 29, 30], the design of adaptive filters [31] 
and the analysis and computation of local and running cross-correlation functions 
using Lagurre cross-correlation sequences [32, 31]. The expression (3.55), although 
formidable in the time domain, will be shown to have an exceptionally convenient 
interpretation in the transform domain which not only provides insight about prop¬ 
erties of the discrete Laguerre functions, but also leads to simple procedures for their 
computation. However, we first show that the Wiener G-functionals can indeed be 
expanded in terms of the discrete Laguerre functions and therefore that the structure 
shown in Fig. 3.1 is an appropriate one. 

2. Orthogonal Expansion of the G-Punctionals 

Given a set of orthogonal functions in one variable, a function of several 
variables can be expanded as a sum of products of these functions. Thus the Wiener 
kernel can be expanded in the p-dimensional series 

oo oo 

^2'* * ' ^ ^p) ^ ^ ^ *** ^ ^ rn2 ••♦TTip ^mi (( ^^2 ) * ' * ) (3.56) 

mi=0 Tnp=0 
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where Xmik) is the degree discrete Laguerre function and the coefficients Crojmj-m, 
are given by 

Cmimj-tn, ~ ' 5^ 9p{^U ^2? ‘ ? kp)Kni {^'l)^Tnai^2) ' ‘' (3.57) 

fc,=0 *,=0 

In this expansion the kernels remain symmetric with respect to the time arguments 
since the indices mi through rrip aissume the same set of values from the the integers 
in the set [0, oo). By Applying the linearity property (3.22) and the expansion in 
(3.56), it can be seen that the p*'* degree G-functional has an expansion of the form 

OO OO 

Gp[gp-,w{n)]= ••• 5] Cm,m 2 -m,Gl[Xm,Xma---Xm,;w{n)] (3.58) 

mi=0 mp=0 

where Gl is a nonhomogeneous functional with leading kernel 9p — <^7111 ' ■ ‘ Xfn ^. 

The kernels of each term in this functional expansion may have the Laguerre functions 
repeated; therefore we write an individual kernel in the form 

flt = \pi \P2 ...XPN (3.59) 


such that 


J2Pj=P A' = 1,2. •••or p 

i=i 


(3.60) 


It is important to point out that this kernel is not symmetric with respect to its 
arguments; but the sum (3.56) or (3.58) produces a kernel Pp which is symmetric. 
This happens because other kernels in the sum have time arguments that permute 
in the functions To give an example, for N = 2 and pi + P 2 = p the leading 

homogeneous functional Gp[pp:u’(n)] can be written as 
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J ... f; 

fci=o k,=0 

xw{n — ki)w{n — ^2) • • • w{n — fcp) 

*1=0 *pi=0 

= G;,lASf'>;<»(n)lG;,lAS^>;u)(n)l 

To keep the symmetry property we could repeat the expansion terms (pi d-Ps)! times 
permuting the time arguments and divide their sum by (p. + Ps)!. In this case 
however, each term would give us the same value, namely 

G„lAif‘>;u.(n)lG„|Ai”';u-(u)l 

Then we would multiply again by (p. + Ps)! which leads to the same result if we 
assume that the kernels of the homogeneous functional are symmetric. In general, 


for any value of A 


G*lAlf>Alr>---AS'''>;u'(u: 


= G;,l,\Sr>;tc(n)lG;.[A!”>;u-(nj 


•Gj„[AS'''>lu.(n) 


(3.62) 


The derived kernels of the G-functional are related to the leading kernel by (3.14)- 
(3.17) and the functions \i are orthogonal i.e. 

f; \i{k)\i(k] = iii (3.63) 

fc=0 

In the above example, the derived homogeneous functionals are obtained by reducing 
the order of the kernel using (3.17). letting t^vo time arguments be equal and summing 
over their value in the range [0. oc). Since the symmetry of the kernel includes terms 
of all possible permutations of the time arguments, we have three cases 



• Summing the product Xi^{k)X^^{k) which equals one; therefore the exponent pi 
is reduced to pi — 2 while p 2 remains the same. 

• Summing the product Xij{k)Xij{k) which equals one; therefore the exponent pj 
is reduced to p 2 — 2 while pi remains the same. 

• Summing the product Xi^{k)Xi^{k) which equals zero from the orthogonality 
property. 

Thus the required functionals are obtained by just reducing the exponents two at a 
time and permuting all of these reductions among the A s. In this way we arrive at the 
results mentioned in (3.26) for the G-functional in which Qp[gp: u'(n)] is represented 
by sums of products of terms t/i(n) as given by (3.27). 

3. Realization of Section 1 of the Wiener Model 


Having shown that the Wiener G-functionals can indeed be expanded in 
terms of the discrete Laguerre functions, let us return to the form of these functions 
and the corresponding realization of section 1 in Fig. 3.1. Equation (3.52) can be 


put in the form 




p[~- p M 


or 


A^-) = Ao(-^) 


(-'-/>) [ (~—/^) J 

P[z- P~'^) 


I iz-p) J 


for which 


Ao(r) = 


(3.64) 


(3.65) 


(3.66) 


{z-p) 

Since for |p| < 1 the term in brackets is the transfer function of an allpass causal 
stable linear filter 


Hap{ = ) = 


P{z-P^) P-Z 


[z-p] 


(3.67; 


72 




Figure 3.4: Realization of section 1 of the model, 
it follows that Ai( 2 ) can be written as 

Ai{z) = A,_i( 2 )i/„p( 2 ) (3.68) 

or 

Ai(.-) = Ao(-')W:,(2) (3.69) 

Since section 1 of the model required a bank of the Xi(n) this means that section 1 
can be realized by a linear stable causal system with a single pole at p followed by a 
train of cascaded allpass filters each with a transfer function Hap{z) as shown in Fig. 
3.4. The frequency response of the leading Laguerre system Ao(e^") is shown in Fig. 
3.5 for p — 0.5. while the frequency response of the allpass filter is shown 

in Fig. 3.6. The impulse response of the allpass filter is given by 

hapi'n) = P^{n) — (1 — — 1) (3.70) 

Beginning with 

Ao(A-) = - pyu(A-) (3.71) 

the remaining functions can be computed by the convolution form 


Aj(n) = Ai_i(n) * hapin) 


(3.72) 
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Figure 3.5: The leading Laguerre function Ao(e^") for p = 0.5. 



w rad/sec 



Figure 3.6: The allpass filter transfer function for p = 0.5. 
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Figure 3.7: Example of the discrete Laguerre functions for p = 0..5. 
or by the recurrence relation 

Ai(n) = p\i{n - 1) + p\i-i{n) - - 1) (3.73) 

An example of four discrete Laguerre functions for p = 0.5 is shown in Fig. 3.7. Note 
how the energy in the sequence tends to move away further from the origin as the 
order increases. Thus when these functions are convolved with an input waveform 
the higher order Laguerre filters tend to bring in more of the “past history” of the 
waveform. 

D. CROSS-CORRELATION FUNCTIONS AND CROSS¬ 
SPECTRA OF THE OUTPUTS OF THE LAGUERRE 
SYSTEMS 

An important part of this dissertation is to derive expressions for the higher 
order statistics (moments, cumulants, polyspectra) of the output x{n) of the Wiener 
model. Given the structure of the model, these higher order statistics will involve 
cross-moments between the outputs of the linear systems comprising section 1. Since 
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these outputs are Gaussian and have zero mean, it is sufficient to describe their 
cross-correlation. It will be seen that because of the considerable structure brought 
to the problem by the use of the discrete Laguerre functions, these cross-correlation 
functions inherit a number of interesting properties. This in turn provides a basis for 
the organized computation of the Wiener model’s higher-order statistics. 

1. The Laguerre Cross-Correlation Sequences and Their 
Transforms 

The cross-correlation function between the outputs yi^ and i/i, of the linear 
section is given by 

+0} 

[fci=0 *3=0 

fel =0fcj=0 

since w(n) is white Gaussian with variance al 

£ {u’(7? — ki)w{n + / — A' 2 )} = — k2) (3.75) 

and (3.74) becomes 

= <7^£4(/r)Ai,(/ + ^’) 

fc =0 

= (3-76) 

where 

'■<.«(/)= + (3.77) 

fc =0 

The quantity will be called the Laguerre cross-correlation sequence. From 

well-known properties of the ^-transform it can also be expressed as 

(3.78) 

2-j Jc\ 
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The complex cross-spectral density function ( 2 ^) is the z-transform of the cross- 

correlation function then given by 

where 5ii,ij(2) is the ^-transform of the Laguerre cross-correlation sequence 

(3-80) 

with 

ru.M=7^'f ( 3 - 81 ) 

Zttj JCl 

The quantity will be referred to as the Laguerre cross-spectral density func¬ 

tion. Then by comparing (3.78) and(.3.80) and substituting (3.65) we have 

= Ai,(;::)Ai,(r~^) 


(-1 ^2\ ^ 

'p{= - p ^)' 

*2 .-1 


\p{z-^-p-^)] 


(*--/>) . 

1 

iz-^-p 

’piz-p-^}] 

1 *2 

(--■-/’) J 

!“ti 

p [z-p){z-p-^) 

. iz-p) . 



This function has the interesting property that it depends only on the difference 
between the orders of the two Laguerre linear system outputs and not on the actual 
values of the orders. If the order difference is denoted by 

d = 22 ~ *1 (3.8.3) 


the Laguerre cross-spectral density function is 

(l-p^*) 


Sd(2) — 


p {z-p){z-p-^) 


pi- - p ) 


(z- p) 


(3.84) 


and the complex cross-spectral density function of the two outputs is given by 


By the same considerations, the Laguerre cross-correlation sequence is a 
function only of the difference d = t 2 - *1 and can be denoted by 


n.M - ri,.M = uil) 


(3.86) 


Therefore (3.76) can also be written as 

=<^o(0 

Since the Laguerre cross-correlation sequences and their transforms depend only on 
the order difference d, it is interesting to examine their behavior for different values 
of this parameter. The leading cross-spectral density function for which d = 0 is the 
complex power spectral density function of the output of any of the Laguerre systems 


and is given by 


Sy{z) = <7oSo(-) = 




The corresponding autocorrelation function is 


(3.88) 


Ry{l) = (rlroil) = <7^/1 


(3.89) 


Further (3.84) indicates that 

= so{z)Ht,{z) (3.90) 

Fig. 3.8 shows a plot of the sequence ro(/) for p = 0.5 and the corresponding power 
spectral density function 5o(e^'^)i which is 5o(2) evaluated on the unit circle. Starting 
with ro(/) the Laguerre cross-correlation sequence can be derived by the convolution 

form 

rj(/) = rj.,(/).A.p(/) (3.91) 

or by the recurrence relation 

ra(l) = prjiil - 1) + prd-iU) - rd-i{t - 1) (3.92) 
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Figure 3.8: The Laguerre autocorrelation sequence ro(/) and power spectral function 
so(e^^) for p = 0.5. 



Figure 3.9: The generation of the Laguerre cross-correlation function. 


Both of these relations follow from considerations similar to those in section C. Fur¬ 
ther, the sequence can be thought of as the “impulse response” of the “system” shown 
in Fig. 3.9. This picture does not represent a real system, but only a convenient way 
to show the relation between the Laguerre cross-correlstion sequences. 

By substituting (3.84) in (3.81). rj(/) has the representation 


rj(/) 





(3.93) 
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with the following characteristics: 


1. For d > 0 and / = 0 


MD) = -(!-. 


(3.94) 


(This is a restatement of the fact that the j/i(n) are orthogonal.) 

2. For d > 0 and / < 0 a change of variables z = v~^ can be made in (3.93). Then 


r.(/) = -d - = » 


(3.95) 


since the integrand has no poles inside the unit circle. In other words, for d > 0 
the rj{l) are strictly right-sided. 

3. For d > 0 and I > 0 (3.93) can be evaluated as 


rail) = -ii-p)P 




d! dz'i' 


j (-l)”»(d+/-r77)! . 

= i'^ - P )—J~ (m. - l)l(d - JT>)\il - m)\^ 


(3.96) 


which provides an explicit expression for rail). 


We also note here that negative values of d need not be considered explicitl}' since 
from the definition (3.77) with d = — ^2 


.Jl) = rA-l) 


Characteristics 1 and 2 above imply that the cross-correlation function (not the auto¬ 
correlation) of the outputs of the Laguerre systems is zero if the output of a higher 
Laguerre order system is averaged with the output of a lower Laguerre order system 
at an equal or greater time lag. In other words, the output of a Laguerre system is 
orthogonal to (uncorrelated with) the future values of the output of the lower order 
Laguerre system. An example of the first four Laguerre cross-correlation sequences 
is shown in Fig. 3.10. 
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Figure 3.10: Laguerre cross-correlation functions for p = 0.5. 

2. The Effect of the Parameter p on the System 
Characteristic 

The above derivations of the transfer functions and cross-spectrum (impulse 
responses and cross-correlation) were carried out for real values of the system pole 
location p. The impulse response of any stable causal system can be expanded in 
terms of the complete orthonormal set of discrete Laguerre functions. The expansion 
generally requires an infinite number of these functions. When an impulse response 
is approximated as a finite weighted sum of Laguerre functions however, the error 
in this expansion is affected by the choice of the parameter p (pole location in the 
complex 2 -plane). Stated another way, for a given mean square approximation error, 
the minimum number of terms in the expansion is reduced when p is appropriately 
chosen. For positive real p the linear system has a lowpass characteristic. If the 
system is required to have a highpass characteristic the best choice of p is on the 
real negative axis. For a system with a bandpass characteristic the pole is best 
chosen as a complex value. This leads to a system with a complex output. However. 


'-5 0 5 10 15 
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combining the response of two systems with parameters p and p* {p conjugate) as 
shown in Fig. 3.11 produces areal output process with essentially the same statistical 
characteristics. This follows because conjugating p merely results in conjugation of 
the functionals Qp in the Wiener representation (see (3.20)). To obtain the real- 



Figure 3.11: The nonlinear model for bandpass processes 

valued output shown in Fig. 3.11 we can equivalently work with the top system and 
take the real part of its output. When a complex value of p is used, (3.66),(3.67) and 
(3.84) need to be modified slightly. 

1. The leading Laguerre function becomes: 

Ac.(.-) = i/l - (3-97) 


2. The allpass filter transfer function becomes: 

w /-A 

{z- p) I- pz ^ 

3. The leading Laguerre cross-spectrum becomes: 




[z - p)(r-i - p*) 


(3.98) 


(3.99) 
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These have the frequency characteristics similar to those for real positive p except 
for a shift on the frequency axis equal to Ip. 
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IV. CUMULANTS AND POLYSPECTRA OF 
THE WIENER MODEL OUTPUT 


As has been mentioned earlier, our interest in the Wiener model is to represent 
a class of discrete random processes. The output of the Wiener model has been 
shown to be represented as the weighted sum of multinomials (3.34) known as the Q- 
polynomials. These multinomials are denoted by Qa^ and each was shown in (3.33) 
to be formed as a specific product of the outputs of the memoryless nonlinear blocks 
of section 3 of the model. Since each of these outputs is represented by a Hermite 
polynomial formed from the output of a linear system driven by a zero-mean white 
Gaussian process, the Q-polynomial is a sum of functionals of a white Gaussian 
process similar to those discussed in Chapter II. Thus the procedures developed 
in that chapter to calculate the cross-moments and the cross-spectral functions for 
this type of functionals can be applied here to develop cumulants and polyspectra of 
processes represented by the Wiener model. Moreover, since the impulse responses of 
the linear filters in section 1 are the discrete Laguerre functions, the cross-correlation 
functions and cross-spectral functions of the outputs of the linear filters were shown 
to have special properties which result in a structured method of representing the 
required cumulants and polyspectra of the modeled process. 

In the following section we present a detailed view of the Wiener model and 
define the Q-polynomials for a model of specific dimensions. In the next sections we 
formulate the process cumulants and develop an organized procedure for their eval¬ 
uation. In the last section that we formulate expressions for the process polyspectra 
and develop a procedure for their evaluation. 


So 



A. THE DETAILED WIENER MODEL 


We consider here a finite version of the Wiener model for discrete nonlinear 
systems presented in Chapter III defined by the two model dimensions Nl and N/f- 
The integer Ni is the highest order of discrete Laguerre functions representing the 
impulse responses of the linear filters comprising section 1 of the model; that is, 
section 1 consists of + 1 filters with impulse responses given by the discrete 
Laguerre functions Ao(rj), Ai(n), • • •, The integer Ns is the highest degree of 

nonlinearity in each memoryless nonlinear block comprising section 2 of the model. 
Ns is therefore the order of the highest order Hermite polynomial in these blocks, 
which is the highest degree of nonlinearity in this representation. The complete 
structure of this model is shown in Fig. 4.1 for Ni = Ns — 2. In general if a 
system is used to model a process with nonlinearity degree p. the lower degrees of 
nonlinearity are also included in the representation. This follows since the system is 
represented using the Wiener G-functionals Go-Gi-"' -Gp-i-Gp- The following rules 
are therefore necessary to define the Q-polynomials. 


1. The kernels of the G-functionals are symmetric with respect to their time ar¬ 
guments. Therefore the expansion of the leading kernel in terms of a product 
of discrete Laguerre functions must include all the possible permutations of 
products of these functions with total sum of exponents equal to the order of 
the kernel. This means that the G-functional is a sum of terms each formed 
by one of the possible permutations of products of Hermite polynomials whose 
sum of orders is equal to the functional order. Therefore in the model represen¬ 
tation the p*^ degree nonlinearity is represented by all the possible permutations 
of products of the outputs of section 2 that result in multinomials of degree p. 
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section 1 section 2 section 3 


) 


Figure 4.1: The detailed Wiener model for Ax = 2 and Nn = 2. 

2. The relation between the derived kernels and the leading kernel of the G- 
functional implies that when a Hermite polynomial of one of the linear section 
outputs is included in the product, no other polynomials of the same output 
are included. Therefore in each permutation of the products of the outputs of 
section 2, one and only one output is taken from each nonlinear block at a time. 

3. If the nonlinearity of the process is required to be of degree Ns- the expression 
for the nonlinear system output does not include functionals of order greater 
than Ns- This means that the sum of the orders of Hermite polynomials in 
each term does not exceed Ajy. 















Recall from Chapter III that the Q-polynomials, which replace the G-functional 
representation of the output of the model, are formed as follows: 

1. Each polynomial is assigned a vector of indices 

oci = 

The value of each component Qy € [0,1,2, • • •, Ns] indicates the order of the 
Hermite polynomial of the corresponding output yj{n) that is included in the 
product, i.e., 

Qq!,(^) = QaioQir"OiNi(^) 

= Ha<o(yo(n))Ha„(j/i(n))---Ha.j^^(yAri,(n)) (4.1) 

2. The set of Q-polynomials are weighted and summed over all the possible values 
of the vectors ctj. To keep the symmetry property of the functional kernels, 
each aij assumes the values 0,1,2, • • •, Ns for same vector oci in the set. The 
maximum degree of nonlinearity is maintained by keeping the sum of the com¬ 
ponents of oci less than or equal to Ns- Hence the Q-polynomial representation 
of the output for specific dimensions Ni and Ns is 

= £ caiQaiin) (4.2) 

E[ai]=o 

The quantity E[ai] (see (2.90)) is the order of the functional represented by the 
Q-polynomial Qai(”). When Il[ai] = 0 we have the multinomial Qo which is the 
product of all lio{yi{n)). Since all Ho(j/i(n)) are equal to 1 it follows that Qo is also 
1. This term represents the zero-order homogeneous functional, which is a constant 
value. The weighted sum of Q-polynomials for which E[ai] = 1 represents the linear 
portion of the representation and is equal to the non-homogeneous functional of 
first order Qi[w{n)]. Similarly the Q-polynomials for which Ii[aj] = 2 represent 
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the quadratic part ^2[«^(«)]5 and so on. In our development of the cumulant and 
polyspectra expressions we will need to define a procedure for ordering the indices 
of the polynomials. We choose the following rules of ordering because they result in 
an organized structure of the arrays of the values of cross-moments and cross-spectra 
of the Q-polynomials required in our formulation. The indices are ordered using the 
vector notation 

Oi > OLj (4.3) 

with the following interpretation: 

1. Oj > oij if S[a,] > S[Q!,]. In other words, the former is greater than the latter 
if it represents a functional of order greater than that represented by the latter. 

2. If Sfoii] = SfcTj], then oii > ocj if the partial sums (below) satisfy 

Nl Nl 

S 3m € [0,1,2, (4.4) 

lc=m k=m 

In fact (4.4) defines the entire ordering since it includes the case m = 0. This 
means that if the Laguerre systems are divided into two groups, one vector of indices 
is considered greater than the other one if the group of the higher order Laguerre 
systems is represented by Hermite polynomials that have a greater sum of orders. 
For example if Nl = 2 and Ns = 2 the output x{n) is comprised of the weighted 
sum of ten Q-polynomials that have indices with the following descending order: 

002 > Oil > 020 > 101 > 110 > 200 second order (quadratic) terms 
001 > 010 > 100 first order (linear) terms 
000 zeroth order (constant) term 

(4.5) 
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Nl/Nk 

1 

2 

3 

4 

1 

3 

6 

10 

15 

2 

4 

10 

20 

35 

3 

5 

15 

35 

70 

4 

6 

21 

56 

126 

5 

7 

28 

84 

210 


TABLE 4.1: Number of section 3 parameters Coj for a set of model dimensions. 

In this example x(n) can be represented as 

002 

H caiQa,- (4-6) 

ai=ooo 

where the coefficients cot- are the parameters that define section 3 of the model. The 
number of coefficients A' depends only on the value of the model dimensions Ni and 
Nn and increases with both dimensions, as shown in Table 4.1. 

B. CUMULANTS OF THE OUTPUT OF THE WIENER 
MODEL 

If the output of the Wiener model is represented as a sum of Q-polynomials (see 
(3.34)) then the order output cumulant can be written as 

' ik-l) - 

cti aj a* 

xcum(Qaj(»i),Qa3(n + h),-” ‘Qoiki'^ + k-i)) 

(4.7) 

To be notationally correct we should denote the a’s by but this nota¬ 

tion gets extremely cumbersome in the development that follows. Therefore for 
example should be interpreted as a generic index vector not simply the first one in 
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the set, and the sum over is the sum over all such index vectors. We can write 
(4.7) more simply as 

Ci'-\l) = Y,caCff(l) (4.8) 

o 

where 

= cum(Qai(n),Qaa(« + 
and a is a concatination of the vectors Oj, that is 

Qt = [a/, aj, • • •, txk^f (4.9) 

while ca is the product of the coefficients 

ca = ca^ca,-”cak 

and the sum in (4.8) is over all the possible permutation of the indices. 

The number of terms in (4.7) would seem to be equal to the number of the 
coefficients cq raised to the power k. which is a huge number even for low dimen¬ 
sional models, and increases exponentially as the order of the cumulant increases. 
Fortunately this is not the case because the cumulant C^\l) of Q-polynomials has 
properties that make its value identically zero for many well-defined combinations 
of the indices. In the following, we define the combinations of indices for which the 
terms in (4.7) are not equal to zero. Moreover from the combination of indices we 
develop a formula to compute the value of the cumulant. 

To obtain an expression for the A-'^-order output cumulant of C^\l) we need to 
develop an expression for the A‘^-order cross-cumulant of the Q-polynomials Cq\1). 
For these cumulants we need to compute the expectation of the product of the Q- 
polynomials (moments) and substitute in the moment-cumulant relations described 
in Chapter II. To develop the expression for these expectations we expand the Q- 
polynomials as the product of Hermite polynomials. These are in turn functions of 
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the Gaussian outputs yi of the linear filters. The values of the components of the 
vectors of indices a together with the time lag arguments of the Q-polynomials can 
be used to determine if in a product of Laguerre cross-correlation sequences there 
exists at least one which is zero. This principle is used to construct a set of sufficient 
conditions on the vectors of indices such that the cumulant of the Q-polynomials 
is not identically zero. Then we can apply the procedure developed in Chapter III 
to find the desired expression. The analysis is presented in three subsections. First 
we form a general expression for the moments of Q-polynomials. Next we show by 
example how we simplify the general expression and eliminate most of the terms. 
Finally we present the general efficient procedure for computation of the cumulants. 

1. General Expression for Moments of Q-Polynomials 

For a model of dimensions Ax, and AV the cross-moment of the Q-polynomials 
can be expressed as 


== £^{Qaa(«)Qa,(« +/i) •••Qafc(» + 4-i)} 




E[m]^! /^2\ 

-2m)! V2 j 


Sfm] 




(4.11) 


where 7?aUm(^) is defined as the moment 

+ h) ■ •. Y“*-'“‘(n + 4_,)} (4.12) 

and where the following vector notation is used 
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m = 

T 

mi = [rriio, mu, • • •, 

i?j l^J iTi 

E = E E ••• E 

xn=0 moo=0»»oi==0 0 

Y“i-’”‘(n) = !/f"'’^(n)!(r'“’”“(")•■•!'«"'■ * 

We seek here to simplify this general expression, to outline its specific structure, and 
to show which terms need to be computed and which do not. The cross-moment of 
the Q-polynomials has been expanded as a weighted sum of expectations of prod¬ 
ucts of Gaussian variables vg’- The weights in this summation are products of the 
coefficients of the Hermite polynomials. Since the sum of the exponents in each 

cross-moment 77 a- 2 m is 

vr« _ 9ml = Efal - 2S[ml 


then the sum of exponents is even if Sfa) is even. Thus a necessary condition for the 
cross-moments of the Q-polynomials to be non-zero is that the total sum of indices 
of the polynomials Ela) is even. The value of the cumulant in (4.8) is therefore the 
sum of l/ir terms uM indices a for o-hich each -(oil is even. In this case the value 
of the cross-moment is obtained by following the procedure described in Chapter 
II. The symmetry properties of the cumulant functions of teal stationary random 
processes are utilized in this respect to calculate the cumulants in the non-redundant 


region defined b\’ 


0 > /i > • • • > h-i 


(4.14) 


This choice ensures that the time lag arguments of the submatrices of Cf”' defined 
as R(u) in (2.88) are non-negative. The entries of these matrices ate the cross¬ 
correlation functions of the outputs of the Laguerre filters comprising section 1 of 
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the model and this gives the matrices specific important structure. Recall that the 
Laguerre cross-correlation functions have the following properties: 

1. The index d of each cross-corelation function is a single integer representing 
the difference between the Laguerre orders i and j of the correlated outputs 
{d = j-i). 


2. The cross-correlation functions (not the autocorrelation functions) are strictly 
single sided; for positive values of d the function rj(/) is non-zero only for 
strictly positive values of the time lag /. 

It follows from the above that the diagonal (Ni -I-1) x {Nl + 1) blocks denoted R(0) 


in (2.87) are the identity matrices because for / = 0 

, , / 1 for d = 0 

'"<i(0) = I Q for d 0 


(4.15) 


For cumulant values in the non-redundant region, the off-diagonal blocks above the 
main diagonal are upper triangular and Toeplitz because for d > 0 and / > 0 rd( —/) 
is identically zero and rj(/) in general is not zero. The blocks have the form 

'roil) r^il) r2il) ••• r^,(/) 

0 ro(/) ri(l) rjv^_i(/) 

R(/) = 0 0 ro{l) ••• rjVi-2(/) ; />0 (4.16) 


0 0 


Since the overall matrix 's symmetric the off diagonal blocks below the main 

diagonal are lower triangular. Then for the case of the k^^-order cumulant of the 
output of the Wiener model the matrix of cross-correlation sequences becomes 


c,(/) = 

R(0) R(/i) R(/2) ••• R(/fc-i) 

R^(/i) R(0) R(/2 -/i) ••• R(/fc-i-/i) 

R^(/fc-i-/i) R^ik-i-h) R^ik-i-h) ••• R(0) 

(4.17) 
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where the R(0) blocks, as previously noted, are the identity matrix. 

The corresponding + 1) x k{NL + 1) multiplicity matrix M is con¬ 

structed such that no permutation resulting in a zero-valued term is allowed. Since 
the moment expression is a sum of the products of system cross-correlation sequences, 
the value of any term including a zero-valued entry of is zero. In any permutation 
we must keep the multiplicity of a zero-valued entry of Cy equal to zero. Then an 
entry of M is zero if the corresponding entry of Cy is zero. This means that the 
multiplicity matrix M has the same structure as Cy, i.e., 

1. The matrix M is partitioned into k x k square blocks each of size (Al -f 1) x 
(A^l + 1)- 

2. The diagonal blocks are also diagonal containing non-negaUve even-valued en¬ 
tries along the diagonal. 

3. The upper non-diagonal blocks are upper triangular and the matrix is symmet¬ 
ric. 

4. The values of the entries of M in the last column and in the last row or along the 
diagonal of the lower right block are dependent on the values of the preceding 
entries in the columns or rows. This dependence is required to maintain the 
condition that the sum along any row or column is equal to the corresponding 
exponent - 2rnij. V^’e write this more simply using the reduction notation 

as 

T[M] = a - 2m 

As an example, consider computing the second order cumulant of a model 
with dimension Nl = 2 • The matrix of the Laguerre cross-correlation sequences has 
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the form 


1 0 

0 1 

0 0 

ro(0 0 


(4.18) 


roil) 0 0 1 0 0 

ri(/) ro(/) 0010 

.^2(0 ra(/) ro(/) 0 0 1. 

where the above mentioned properties are evident. The associated multiplicity matrix 
M has a similar form and must satisfy the condition 


roil) rril) r,il) 

0 ro(/) ri(/) 

0 0 roil) 

1 0 0 

0 1 0 

0 0 1 


2ii 0 0 is 14 

0 2is 0 0 io if 

0 0 2^8 0 0 *9 

h 0 0 2?:io 0 0 

?3 Js 0 0 2?ii 0 

*4 i? *9 0 0 2ii 


QTlO — 21X1x0 
axi - 2mix 
Qi 2 - 2mi2 
020 ~ 27X120 
021 - 2m2i 
O 22 ~ 2x1122 


(4.19) 


Note that the vector on the right hand side represents the exponents of yi in a 
specific term in the expression (4.11), and that the moment 7?2Um(0 is expressed 
using (2.92) as 


^a-2m(0 = 


siCKj-jsfin] 


(<^o) ’ (a^ - 2 mi)!(a, - 2m2)! 


V 1 

2 [M]=a — am 2 » \ 2 )■ h=3i+i 


MLhJil 2 ( 2 + 1 ) 


iCyijuj2))^’'^^'^'> 

MijxJfV. 


(4.20) 


In this example the entries of M are divided into two groups; the first is that of 
the independent entries {ii,is-is-is-is-is} and the second is that of the dependent 
entries {i4, i?? isi iio-in? iis}- The summation 

E[M] = a - 2m 

is actually a multiple summation over the non-negative integer values of the inde¬ 
pendent entries while the dependent ones are determined from condition (4.19). In 
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general, for a system of any size Ni, the value of the fc*^-order cross-moment is given 


rrCK-am 


»?S- 2 m (0 = - 2ma)!(a, - 2mj)! • • • (a* - 2m*)! 

^ 1 *<%«> (C,0..i,))“^ “"A*” (C,(j., 


(4.21) 


Now we are ready to consider the cross-moment of the Q-polynomials which 
we denote by A necessary condition for this statistic to be non-zero is that S[q:] 
is an even value and the permutations are arranged so that the specific structure of 
the multiplicity matrix M is kept as described above. In this case the value of the 
cros-moment is calculated by substituting (4.21) in (4.11) to obtain 




L^J 


(-l)^t“laJaJ-.-afc! 


~ n^o 1^0 nii!m 2 ! • • • mfc!(ai - 2mi )!(a, - 2 m 2 )! • • • (a* - 2m*)! V 2 


x((7^)*’°'’»**”'' (Oi - 2mi)!(aa - 2 m 2 )! • • • (a* - 2m*)! 


XI .,(M) n 

E[M]=Q'-2in 2 2 ii=l 


n 

h =ii+l 


A/(ji,j2)! 


Notice that since the first summations are over the terms mi, the dependent entries 
of the multiplicity matrix change to maintain the condition S[M] = a — 2m. The 
last equation can be simplified to 




l^J 1^1 L^J 


V- 1 ))««•«) 


^ s .,iMi n 

E[M]=:Q!-2ra 2 2 Ji=l 


Ji ^3l) \l 
2 


i2=ji+i 


Af(ji.;2)! 


It appears that the general expression (4.23) for the cross-moment of Q- 
polynomials is extremely complex. However it will be shown that considerable sim¬ 
plifications are possible and that when moments are combined to produce cumulants, 
even further simplifications arise. These combine to make the overall computation 
quite managable. 

2. Rules for Computation of Moments of Q-polynomials 

The key to reducing (4.23) to a simple form is the following. Although we 
have specified some necessary conditions for the terms in (4.23) to be non-zero, these 
conditions are not sufficient. In fact, if we proceed in calculating the cumulant of 
the Q-polynomial according to the arrangement in (4.23) many computations are not 
necessary due to the following: 

1. For a large number of configurations of the values of the first summation pa¬ 
rameters m,- there is no multiplicity matrix M with non-negative integer entries 
that satisfies the the condition S[M] = a — 2m. In these cases the correspond¬ 
ing term is equal to zero. 

2. In many other configurations there exsits a matrix M satisfying the above 
condition but the resulting terms sum to zero. Therefore the value of the 
moment is also zero. 


This results in a tremendous simplification of the resulting moment expressions and 
a corresponding savings in computation. For instance, in the case of the second-order 
moment of Q-polynomials of our example, (4.23) reduces to the considerably simpler 
expression 




= {al) » aja,! ^ 

» 3=0 


aio-a2o **2+»«+»9/ 


(/)ri*+*^(/)r‘M/) 


f2fi3h4!?7!?9! 


(4.24) 
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where the quantities and *9 depend upon the value of *3 as described in 

Appendix C. 

In Appendix C we examine the conditions for which we can avoid computing 
the terms in the Q-polynomial given by (4.11) which are identically zero. This leads 
to the following principles; 

1 . The non-zero terms in the expression for the cross-moment of Q-polynomials 
come only from the expectation of the product of the first terms in the Hermite 
polynomials corresponding to m = 0. This means that the Q-polynomial cross¬ 
moment expression ( 4 . 11 ) reduces to a single term and that the computational 
effort is correspondingly reduced by a large factor. 

2. This expectation is taken in a special way such that the pairing permutations 
do not include any permutation that would result in the term ro(0) as a factor. 
This corresponds to setting the values of the diagonal entries of the multiplicity 
matrix M equal to zero. 

Although the Q-polynomial expression (4.2.3) is reduced to one term only, this term 
is also investigated to determine the relations between the components of the index 
vectors cti such that this term itself is not identically zero. The relations between 
the index vector components depend upon the order of the moment: therefore we 
develop these relations for each moment order below. 

a. The First Order Moment (Mean) of the Q-polynomials 

The procedure described above can be applied to compute the value of 
the first order moment. In this case, according to the principles just cited, we take 
the moment of the product of only the first terms of the Hermite polynomials without 
allowing the pairing that results in having ro(0) as one of the quantities within the 
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product. In this case however, the matrix of cross-correlation sequences Cy in (4.17) 
is the upper-left block denoted by R(0) for which all of the off-diagonal entries are 
zero and the diagonal entries are ro(0) = 1. Also the M matrix is diagonal with 
even valued diagonal elements. If the sum of the components of the vector a is odd 
then it is clear that the mean is zero. For any a with even sum of components, the 
diagonal elements of M must satisfy 

E[M] = a 

This means that if a / 0 the diagonal elements are not equal to zero and each 
of the individual components Qjj of the vector tx must be even. In other words, 
when expanding the Q-polynomial using the Hermite polynomials and taking the 
expectation, the sum of the resulting terms has the form 

(1 - 1)^(1 - - 1) » (4.25) 

which equals zero except for a = 0 . In that case 

,4'> = £Wo(n)}= 1 

Let us now apply this result to compute the mean of the process x{n) which is given 

by 

£{x(n)} = £■ | 5 ^caiQai(») 

= (4.26) 

oci 

All the terms under summation are zero except the first one for which 

S {^(i?)} = CoS {Qo} = Co (4-27) 

Therefore cq is equal to the mean of the process. Since the cumulants of a random 
process do not change when a constant is added to or subtracted from the process 
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we will henceforth assume that c© is equal to zero and thus that x{n) has mean zero. 
If Co is not zero, this affects only the mean of the process and not the higher order 
cumulants. 



where the operator Sh denotes the expectation obtained by summing up all the 
possible pairing permutations excluding those that result in one or more instances 
ro(0). Since 

<T^ro(0) = 5{yi(n)yi(n)} = S{yi{n -f /)y,-(n + /)} 

this means that all of the yj{n) or y,-(n + /) of the first term of a Hermite polynomial 
must be paired with those from other polynomials. Note also that the pairing with 




other yj that have the same time argument is equal to zero because of the property 
rrf{0) = 0 for c? ^ 0. Finally, note that if is / is positive the variable yi{n) of some 
Laguerre order can only be paired with another variable yj{n) with Laguerre order 
j > i. This follows because rrf(/) = = 0 for j < i. Therefore, let us arrange 

the components of the index vectors according to the Laguerre order and the time 
index as follows 


yo 

yi 

y2 

••• ytvx ,-2 

yNi,-i 

ysi. 

n + l: O20 

021 

022 

• • • Q;2(JVj.-2) 

a 2 (tVi-l) 

0^2Nl 

n : oio 

On 

012 

• • • Q:i(Ari,- 2 ) 

Ol(^i-l) 

OlJVi 


To obtain a non-zero result yjVi(n) can be paired with -f- 1) only. Therefore the 
exponent of the latter must be greater than or equal to the exponent of the former, 
i.e., 

CtZNj, > CllAfi (4-31) 

Also since ?/;v^_i(n) can be paired only with yp/^-i{n -|- /) or + 1) to obtain a 

nonzero result, the sum of the exponents of the two outputs with highest Laguerre 
orders must satisfy 

Q2Ari + Q2(jVi-i) > QiATx, + Qi(JVx.-i) (4.-32) 

Continuing this argument leads to the set of conditions 

+ Q2(JVt-l) + Q2(JVi,-2) > OilNt + <^l(Ari,-l) + 

+ (>2(Ni,-2) + ^2{N[.-3) > + Ol(Afx,-l) + <^l(JVi-2) + 0'l(;\ri-3) 


O2N1, + C»2(Ni,-l) + " • + 022 + O21 > <^ 1 Nl + Oi(iVi,-l) + • • • + O12 -f On 

O2N1, + 02(JVi-l) + • • • + O22 + 021 ■+ O20 ^ Oi^^ -{■ Qi(;v^_l) -!-••• + O12 -t- On -f Qio 

(4.33) 
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Now let the vector be defined as the vector composed of the rightmost m compo¬ 
nents of the vector a^. 

0i = • • •, oiiNi ,, (4-34) 

This means that /3® represents the exponents of the outputs of the q systems with 
highest Laguerre orders. With this definition, (4.31)-(4.33) can be summarized as: 

{oTq=h2,---,NL + l (4.35) 

where the last relation for q = Ni-\-1, implies that the total sum of the components 
in tta is greater than the total sum of the components in ttj. 

Now. let us consider the relation of the index components starting with 
the outputs of the lowest Laguerre orders. The output ?/o(” + 0 be paired w’ith 
t/o(n) only. Therefore 

oio > 020 (4.36) 

The output yi{n + 1) can be paired with yo(n) or yi{n) only, which implies that 

OlO + Oil ^ 020 + 021 (4.37) 

By similar consideration the following relations must hold: 

OlO + Oil -|- Qi 2 ) > 020 + 021 + O22 

OlO + Oil d" O12 + 0'13 > O20 + O21 -H O22 + O23 


OlO + Oil +-h Oi(iv ^_2 + Oi(;V^_l) > O 20 + O 21 -h-1- 02[Nl-2) + 02(Arx,-l) 

OlO + Oil -h • ♦ • -h Oi(;v^_l) -f > O20 + O21 -f • • • + Q 2 ( 7 \r^_l) -I- Q2N1, 


(4..38) 
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If the vector -y? is defined to be composed of the leftmost q components of oti (repre¬ 
senting the exponents of the outputs of the q systems with lowest Laguerre orders), 

7? = «fi, • • * 1 (4.39) 

then the relations (4.36)-(4.38) can be summarized as 

7’ ^ 7i for 9 = 1,2, • • • yNi, + 1 (4.40) 

The last relation in (4.40) (for q = Nl + 1) means that E[ai] is greater than S[aa]. 
This and the last relation in (4.35) (for m = Nl+I) are not true unless E[aJ = Sfa,] 

. This implies that: 

1. S[a] is even where a = [aj'. 

2. Since <Xi = [7,^^"*^^”* f^i] for ^riy q the relations (4.40) are true if and only if 
the relations (4..35) are true. 

From this condition we also observe that if the Q-polynomials in the expression (4.2) 
are divided into groups according to the degree of nonlinearity they represent (i.e. 
linear, quadratic, cubic, etc.), the necessary and sufficient conditions for the second- 
order cross-cumulant of tw’o Q-polynomials to be not identically zero are. 

1. The two polynomials belong to the same group of nonlinearity. 

2. The relations (4.33) or (4.35) hold. 

c. Third-Order Moment 

The relations between the vectors of indices a, and oc^ such that 
nS{hJ 2 ) = S{Qa,{^^)Qoc,{n + h)Cla^(n + h)}^0 for/j > /i > 0 (4.41) 
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are developed as follows. The components of each vector are arranged according to 
descending time argument value and increasing Laguerre order of the corresponding 

output: 

yo yi Vi 

n + h- 030 03I O32 • • • OC 3 {Ni,- 2 ) 

n + h: 020 021 022 * • • oc 7 {Si,- 2 ) 

n : oio On ai2 • • • Oi(jv^_2) oi(Ari,-i) oiiv^. 

We have two conditions on these variables. The first condition follows from the 
condition on the expectation of the product of Gaussian variables; namely S[a] is 
even. The second condition is obtained if we examine the conditions on the values 
of the individual components in each vector. We find the relations 

for g = 1,2,• • •,Ni + 1 (4-42) 

^9 _j. > -y? for g = 1,2, • * •, A'i + 1 (4-43) 

+/3® >/3’ for q = 1,2 ,---,A'l +1 (4-44) 

where the vectors and 7 ? have the same definition as before. For q = A'l + 1 
equations (4.42), (4.43) and (4.44) yield the following relations: 

S[a,] < I[a,] + SK] 

S[a3]<S[a.] + S[a,] 

S[a,] < S[a3] + Sla,] (4.45) 

If these three relations are true one of the relations in (4.42), (4.43) and (4.44) 
becomes redundant. Therefore necessary and sufficient conditions for the values of 
the components in the three vectors of exponents to ensure that the value of the 
moment is not identically zero are : 

1 . The sum of all the components of the three vectors !][o£] is even. 

2. The sum of the components of any two vectors is greater than or equal to the 
sum of the components in the third. 
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3. The relations (4.42) and (4.43) are satisfied. 

Similar relations can be developed between the vectors of indices of the Q-polynomials 
for computing the cross-moment of any order. The set of indices in the expression of 
the cumulant of the output of the model is generated then examined to exclude the 
terms that do not satisfy these relations. The value of the cumulant is then obtained 
by calculating the rest of the terms because their values are not identically zero. In 
general, the non-zero valued Q-polynomial cross-moment in (4.11) is then given by 


s[M]=:a j 




(4.46) 


where the first summation is over the possible multiplicity matrix permutations. The 
product is over all the nonzero-valued entries ij of the multiplicity matrix and Vj is 
the corresponding Laguerre cross-correlation with time argument Ij. 


3. Computation of Cumulants of Q-Polynomials and 
Model Output 


To obtain an expression for the A'‘^-order cumulant of the output r(n) as 
described by (4.7) we need to develop an expression for the A'‘^-order cross-cumulant 
of the Q-polynomials 

CgVl.(!•••■, 4 - 1 ) 

Since the Q-polynomials are all zero-mean except Qo (which has a mean equal to one) 
the mean of the model output was shown to be given by the value of the coefficient cq 
(see (4.27)). In our development of the model output cumulants we assume that this 
coefficient equals zero so the model output henceforth has zero mean. Consequently 
the Q-polynomial representation of the model output does not include Qao- The 
expressions of the model output cumulants and the Q-polynomial cumulants are 
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given by (2.14). In this case the Q-polynomial cumulants of orders 2 to 4 have the 
expression 

Cg\l) = ^{Qa,(n)Qa,(n + /)} 

= ^{Qa,(n)Qa,(n + /a)Qa3(n + /2)} 

= 5{Qa,(n)Qa,(n + /i)Qa3(n + /2)Qa3(n + /3)} 

{Qai('^)Qaa(n + k)} £ {QQt3(»^ + l2)Qa^in + I3)} 
{Qai(”)Qa3(^ + ^ 2 )} ^ {Qaa(” + h)Qa^{'n + /a)} 
{Qai(n)Qa4(^ + ^ 3 )} ^ {Qaa(^ + ^i)Qa3(” + 4)} 

(4.47) 

It can be seen that the second- and third-order cumulants of the Q-polynomials 
are given directly by their respective moments. As for the fourth-order cumulant, 
we still can compute it directly using the same procedure we have developed to 
compute the moment. Since the last three terms cancel with similar terms that arise 
in the computation of the first term, we need only to compute the moment using 
the procedure developed above and set a condition so that the common terms which 
cancel are not computed. This condition is imposed on the permutations simply by 
setting terms of the multiplicity matrix M to zero. For the fourth-order moment this 
matrix has 4x4 blocks and we do our computations on the upper trianglar parts 
of the six ofF-diagonal blocks. The entries of this matrix take on values as described 
before such that at least three of the six blocks have entries which are not all equal 
to zero. 

The computation of the model output cumulant can be summarized as fol¬ 
lows: 

1 . From the model dimensions h’l and Nn construct a vector of indices a the 
partitions of which are the Q-polynomial indices Oj. 
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2. For computation of a cumulant of order k, test all the possible combinations 

of k indices a = to determine those that satisfy the 

conditions for a nonzero-valued cross-moment. 

3. For each of the combinations that give nonzero-valued moments use (4.46) 
to obtain the value of the cumulant by applying the proper condition on the 
permutations of the matrix of multiplicity. 

4. From the values of the coefficients coj the output cumulant is given by (4.7). 

C. POLYSPECTRA OF THE OUTPUT OF WIENER 
MODEL 

The polyspectra of the output of the Wiener model are the Fourier transforms 
of the corresponding cumulants of this output. Due to the linearity of the Fourier 
transform these polyspectra are the weighted sum of the polyspectra of the individ¬ 
ual Q-polynomials in the expression (4.2) with the same coefficients ca- Although 
it might at first seem that the procedures developed in the previous section could 
be adopted directly to compute the polyspectra, the regions of symmetry of these 
quantities and other considerations make the problem sufficiently complicated that it 
is better to develop the spectral procedures on their own. The resulting expressions 
and procedures for computation are analogous to those for the cumulants but differ 
in important details. 

Since the Q-polynomials are functionals of Gaussian random variables yi(n), the 
expressions derived in Chapter II Section C.2 can be applied to the compute cross¬ 
polyspectra of the Q-polynomials. The second-order system cross-spectral functions, 
which are the entries of the matrix S defined in Chapter II. are shown in Chapter III 
to be the Laguerre cross-spectral functions In the following development we 
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make use of the functions 


/=>-oo 

which are used in Chapter III to derive the Laguerre cross-correlation sequences. 
These functions, when evaluated on the unit circle, are the Laguerre cross-spectral 
functions^ We also make use of the orthogonality of Q-polynomials, which can be 
detected from the vectors of indices ot, to represent the polyspectra as a sum of only 
the terms that are not guaranteed to be identically zero. 

Since the Q-polynomial moments of orders higher than third are not equal to 
their cumulants (see (4.47)), their Fourier transform is different than the polyspectra. 
We shall call the Fourier tranforms of these moments the Q-spectral functions and 
denote the ^‘'‘-order Q-spectral function as S^{u;). The Q-spectral function has the 
form 


5?(u;) = 


L^J 


L^J 




m 


E 

E[M]=a-2in 


9 2 


( KNi^+i) J I \ 


cja-2m/ N 
conv \ / 


(4.49) 


where M is the multiplicity matrix, which depends on the summation parameters 
m through the relation I1[M] = a — 2m. As described in Chapter II, the quantity 
jg formed as the product of the two quantities : 

1. The product of values of each quantity in the block matrix R(0) that is included 
in the pairing permutation. 

^We note as before, that the use of the term «d(w) rather than 8d(e^") is slightly abusive but 
results in some simplification of notation. 
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2. The multiple convolution of the system cross-spectral functions 5 ^( 2 ) that are 
included in the pairing permutation. These were defined in Chapter II as the 
entries of the off-diagonal blocks S of the matrix S given by (2.106). 


Since the matrix M determines the multiplicity of each of the above quantities, the 
entries of this matrix need to be set to zero whenever the inclusion of the correspond¬ 
ing spectral functions results in a term which is identically zero. 

We have noted earlier (see Section B.l) that the diagonal blocks of the multiplic¬ 
ity matrix are diagonal. Since each of the off-diagonal blocks S" of S is Hermi- 
tian symmetric, if the upper triangle has entries that are expressed as Sd{u:) for 
d — 0, !,•••, Ai, then the lower triangle entries are given by a.'). The corre¬ 
sponding 2 -transforms are Sd[z) and s<f( 2 ~^). For purpose of this discussion let us 
denote a block of the multiplicity matrix corresponding to a block 5“” in the spectral 
matrix as M". Now if the pairing permutation results in a combination of entries 
from the both the upper triangle and the lower triangle then at least one of the 
convolutions takes the form 

^<i,(~) * (“) ^d2{t')dv (4.50) 


If we substitue the value of from (3.84) we get 


SdA = )*^d,{'~ 


2-Kj Jci {v — 2 /)“^)'^!+^ { v — 


_/ (r --pr 

i Jci (v — 


dv (4.51) 


which is identically zero because there are no poles inside the unit circle. Therefore 
if the values of the polyspectra are computed in the non-redundant regions for which 


the off-diagonal blocks M" of the multiplicity matrix M must be upper triangular to 
exclude the permutations for which 5<i(2"^) is included. Then the multiplicity matrix 
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in this procedure is exactly the same as that used for computing the moments. Since 
the terms in the lower trianglular portion are never used, the block S" can be replaced 
by one which is upper triangular and Toeplitz: 

So(u>) Si(u)) S2(u») ••• Siv^(u;) 

0 so(u;) si(u;) ••• 

S" = 0 0 so(a;) • • • (4.53) 

0 0 0 ••• so(<^) 

Notice that a single subscript d = j — i is used for the system cross-spectral functions 
rather than the double index used in Chapter II because the Laguerre cross-spectral 
functions depend only on the difference between the two Laguerre system orders i 
and j. 

In the development in Appendix C of the relations between the vectors of indices 
for w’hich the value of the moment is not identically zero the results were based upon 
excluding permutations that contain the terms ro(0). These are the diagonal entries 
of the matrix Cy in (4.17) and result in a common factor that equals zero. Since 
these functions also appear in the polyspectra expressions under similar conditions, 
the same factorization exists with the same zero values. This means that terms of 
the Q-spectral functions are identically zero in exactly the same cases for w'hich the 
moments of the Q-polynomials are identically zero. The relations between the vectors 
of indices developed in the previous section are therefore the same for the Q-spectral 
functions. Thus the expressions (4.49) simplifies in the same way as the expression 
(4.23) did for the time domain computation. 

Let us now consider the remaining part of the problem, namely computing the 
quantities It will be seen that, like in other calculations involving the 

Wiener model, there is abundant structure present so that the form of the result 
can be predicted from an analysis presented here and the actual computation can be 
implemented by a table lookup scheme. 





1. Multiple Convolutions of the Laguerre System Cross- 
Spectral Functions 

The Q-spectral functions given by (4.49) are obtained by performing multi¬ 
ple convolutions on the system cross-spectral functions 5 ^( 2 ) with given multiplicity 
to compute the frequency-dependent quantities Convolving the spectral 

functions by directly applying the convolution formula (4.50) is a tedious operation 
especially when this needs to be repeated many times. Fortunately however, we can 
use the characteristic of the Laguerre cross-spectral functions to develop a procedure 
to obtain the value of a multiple convolution algebraically. If a system cross-spectral 
function has a multiplicity m. this quantity equals the sum of the entries along the 
corresponding diagonals, principal and non-principal, in the block M" of the mul¬ 
tiplicity matrix M. This means that the contribution of this function s<i( 2 ) to the 
quantity is the result of convolution with itself (m — 1) times. The result is 

then convolved with the contributions of the other functions. Let us use the notation 
5 ^( 2 : p) to denote the a complex cross-spectral density function with parameter p and 
define ^ 5 *( 2 ; p) as 


4"‘^(2;p) = Sd{z-. p) * Sd{z: p) * Sd(z: p) * ^ p)^ 

m ““ 1 convolutions 

= (i.U) 

Using this notation we can also define 

4 ”^(-;p) = <5(2) 

s^P{z:p) = 5 ,( 2 :/)) ( 4 . 55 ) 


We consider various special cases below. 
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a. Multiple Self-Convolution of the Complex Spectral Density 
Function so{z]p) 

When the sum of the entries along the principal diagonal of the block 
M" is non-zero then the system complex spectral density function so{z‘,p) has a 
contribution to the required spectral function This contribution is obtained 

by performing a number of successive convoutions. The Laguerre system complex 
spectral density function has the form 

^ (4.56) 


•so(2;^) = -■ 


The first convolution is given by 


P {2- p){^-P 


= ^o(z;p) * so(z;p) 
27rj Jc 


1 r _ V _dr (4 57 ) 

''27rj Jci (u - zp)(v - zp-'^){v - p){v - p-^) 


Since \z\ < |p“^l it follows that \zp\ < 1|. Therefore the integrand has two poles 
inside the unit circle, at locations p and zp, and by evaluating residues the result of 
this convolution is found to be 

( 2 ). . _ (1 - P^) _£_ 

p2 (2 _ p2)(^ _ ^-2) 

= so(^:p=') (4.58) 


It is seen that the two poles of the complex power spectral density function are moved 
from the locations p and p~^ to and p"^ respectively as a result of one convolution. 
One more convolution can be similarly shown to result in the quantity 

= so[z\p)*sl[z-.p) = SQ{z\p)*so[z:p'^) 

p3 

= so(;:p") (4.59) 
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Continuing in this way we can find that 

4”>(2;/.) = io(2;D (4-60) 

Thus the result of m — 1 convolutions on m functions So(2; p) is equal to the same 
function with the parameter p raised to the m power; its poles are moved from p and 
p~^ to /)”* and 


b. Convolutions of up to Three Laguerre Cross-Spectral Density 
Functions sj^{z;p), sj^{z:p), Sd,(z;p) 

For the other Laguerre spectral functions S(f,( 2 ;/)) for which dj > 0 the 
sum of the elements along the corresponding non-principal diagonal of M" determines 
the number of convolutions to be performed. Let us first consider the convolution of 
with itself: 

^d]\=',P) = -‘‘di(z:p) * Sd,{z:p) 

= ^ / (-)dv 

ZTTJ JCI \V/ 

“ ^ ^ 2;riyci ^ ^ 

Since the integrand has (dj -I-1) poles inside the unit circle at location p the value of 
this convolution is given by 


4, (';/>) = “(1-A 




dildV^^ 


zp-l )(«'>+!) 


1 (rfi-l) (ai-mj-fn2) J J 

m?o ma!m 2 !m 3 !m 4 !(l-mi)! 




7711=0 7712=0 7713=0 

--^- 

7714 =(fi “^1 ”"^2 "“^3 


xtt -i- u^P-P I U -^- r^\p-^p) 

{di-l-wt)] {di-l-mz)l 


di -1-7713 




(di -f 7774 )! _1_ 

di! (p_ ;^-l)(d> + l+m4) 
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(n-,^-m,) ( d,-l\( d^ + mA {-ir*{p-Zp)^'^^-'^-”^^'> 

^ «r^o \ ”*3 / \ / {p- 


m4=n-“mi —mj —tns 


(4.62) 


Since mi can take on values of only 0 and 1, the double summation over mi and m 2 
can be simplified. In particular, the term 


is equal to 1, and with a change of variables 

mi + m 2 = k' 


(4.63) 


k' assumes values either in the range [l,di] (for mi = 1) or the range [0,di — 1] (for 
m-i = 0). Therefore (4.62) becomes 


di inin[fc',rfi-1] 

P)(.;p) = -(l-p’)E E (-1)' 

k*=0 m2=k* — l 

(«'»-*') / ^ - 1 

X 5 : (-ir 

m3=0 \ 


di — mj I 


di-i 


>2(<ii-m2)^2 p2^cii->m2 


dl — I \ ( 2 di — k* — m3 \ 2(2di+l-fc'-m3) 


2(r 

^ {Z - ^2j(2(li+l-fc'-m,) 


(4.64) 


Now let us define 


mi=k'-l V ^ / 


(4.65) 


Since m 2 = 0 for k' = 0. this quantity has the value 




(4.66) 
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Similarly since m 2 = <fi — 1, for A:' = dx we have 




(4.67) 


and for F € [l,di — 1] 




= (-1) 


di^k* ( di 




(4.68) 


Notice that if we let k' = 0 and k' = dj in (4.68) we obtain the relations (4.66) and 
(4.67) respectively. Therefore (4.68) is the required expression for 14'^|^^(A’') for all 
k' e [0,di]. We then substitute k = dx - k' and define as 


= (- 1 )* 


(/>■-!)" I- 


dx - k 


(4.69) 


Hence (4.64) becomes 




di minfAjcii-1] 

V' i 


dx-\ 


dx + A’ — 77 ? 3 
dx 


2(iii+l+fc-mi 




(4.70) 


— p2 j(<ii+1+^ — 7713 ) ' ’ ^ 

The expression (4.70) reveals important characteristics of the convolu¬ 
tion of the two Laguerre cross-spectral functions of the same order. This convolution 
can be obtained algebraically as a sum of spectral functions that have one zero at 
the origin and multiple other zeros at c = 1. The poles of these functions are at 
2 = . The number of these functions depends on the orders n of the cross-spectral 

functions being convolved and they differ only in the multiplicities of their poles and 
zeros. W'e will see presently that similar properties hold for multiple convolutions of 
the Laguerre cross-spectral functions of the same order and convolutions of Laguerre 
cross-specrtal functions of different orders. 
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Now let us define the spectral functions embeded in (4.70) that result 
from the first convolution as 

Then (4.70) can be written as 

di iiun[lB,di~l] 

s^^^{z;p) = E E (4.72) 

fc=0 mj=0 

where 

^3) = (-1 )"■« ( V3 ‘ ^ "O “ ■ 

p^^^k^ms) = di + 2 k 

q^^\mz) = di - 1 - m 3 (4.73) 

with m 3 = 0 , 1 , • • •, /b and A- = 0 , 1 ,2 • • •, di. Therefore the pole and zero multiplicities 
have a range of values given by 

p^^^ = di + T di + 2. • • •, di + 2 — q^^^ 

= 0,1,2, •••.di — 1 (4.74) 

The procedure we have employed so far is useful for the other cases to 
be considered in this section and the results will be analogous. In particular, suppose 
it is desired to form the spectral term = sj, (2) * •sjJ;). We can assume with 

no loss of generality that dz < di. In this case we obtain the same results except that 
the limits on m 3 in (4.72) change from m 3 € [0, di - 1] to 0 < m 3 < min[di - 1, d 2 - A:] 
and the ranges of the pole and zero multiplicities in (4.74) change to 

p^^^ = d\ + 1, di + 2. • • •. di + 2 -f q^^'^ 

q^^'^ e [di — d 2 — 1. di — 1] (4.7.5) 
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Figure 4.2; The multiplicities of the poles and zeros in 5 ^,( 5 ;/)) * Sdj{z]p). 

Fig. 4.2 shows the values assumed by and in (4.72) when con¬ 
volving Sd,iz[p) and 5 ^,( 2 ) for dz < d^. From this plot we note that the minimum 
value of the difference is given by 

min(p^^^ - 9^^^) = min(p^^^) - max(g^^^) 

= (di + l)-(di-l) = 2 (4.76) 


Now let us consider what happens when the spectral term Sdi^dj(-) = -%(^) * 
is convolved with yet another cross-spectral function sjjz). We will show that the 
general structure of the problem observed so far is maintained. The convolution of 
these cross-spectral functions can be written as 


(-)*4iU(-) 


minffctdi —l] 


k=0 m3=0 


b^^\k,m3)sd^{z: p) * X^^\z\\ 


where each term under the summation involves a convolution given by 
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2xj Jci V 

1 f ^p(0-g(»)-i^^ _ p-i)«*»-^(t> - z)^^' 


— i 

2lTj JCi 


{v — p)*‘»+^(v — zp~^y‘^ 


(4.78) 


Since (4.76) establishes that > 1, the integrand has at least one zero at the 

origin and da + 1 poles inside the unit circle at p. Thus we can find an expression for 
the convolution (4.78) similar to that of (4.72), namelj’ 

*=” (4.79) 


where the function is defined as 

v(t), ' p’” --(a - P)» 

X (Z,p,g)- 

and where the remaining quantities in (4.79) are given by 

p(2) _ p(i) _|. jg _ _ ^3 


(4.80) 


- niz 

=(-1,™. ( t;; ) ( 

= ”^E( t ) -1)'*-"-’ (^) 

mo £-1-1 ' ' ^ ' 


q^^'> \ ( — 1 + da — /l 

ma / I da - A’ - 7773 


da - 7772 


mj £-1-1 


A' W da - A: 

7772 / \ da — 777 2 
^t-l-t- 7772 -A-j 


The function has one zero at the origin, zeros at p. and poles at p®. The 
limits of summations in (4.79) determine the ranges of and or equivalently 
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Figure 4.3: The range of the functions as a result of Sd^iz', p) * 
the set of spectral functions For each value of and we have 

p(2) _ p(i)^ p(i) ^ 2 ... p(l) _|. _j_ ^(2) _ ^(1) 

= (4.82) 

Fig. 4.3 shows the region of pole and zero multiplicities of the functions that 
result from convolving 5d3(.:) with a single function X^^^ in (4.77). To define all the 
functions X^^^ that result from the convolution we combine (4.74) and (4.82). Then 
the entire set of spectral functions has pole and zero multiplicities in the region 
defined by 


— di -j- 1, di -t 2. dj -j- ds -h 2 -j- q^^^ 

= 0.1.2. •••.</!-1 (4.83) 

This region is shown in Fig. 4.4. 
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Figure 4.4: The range of the functions as a result of Sd^iz) * 

c. Convolution of Any Number of Laguerre Cross-Spectral 
Ftmctions 

To generalize the results of the previous subsections, let us assume that 
we need to perform convolutions of the Laguerre cross-spectral functions of orders 
c?i, ^ 2 . • •' 1 such that 

di ^ <^2 ^ ^ ^ 0 


We can denote this by 




(2) = Sdi( 2 )* 5^3(2) 

---- — - ' 


K convolutions 


(4.84) 


According to the above discussions we could convolve with then con¬ 

volve the resulting spectral function with s<i,(2) and so on until all the required k 
convolutions are performed. In so doing we can benefit from the above relations to 
establish an efficient procedure to obtain the result. We proceed by induction. 
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Let us define the spectral function X^’'^z]p,q) as one of the spectral 
functions that results from performing v convolutions. This spectral function has 
the general form 


- 1 )’ 
- 1 )P 


(4.85) 


It is seen in the previous subsections that when a function of this form is con¬ 
volved with one of the Laguerre cross-spectral density functions, a weighted sum 
of spectral functions of the same form results. Thus the convolution of Sd{z\p) with 
X^''\z\p^''\q^'''^) produces an expression that can be written as 

Si{z-,p)*X^''\z-,p^''\.q^''^) 

d inin[g^‘'\d—At] . 

= S E h''''^^\k,mz)X^''^^\z\p^^^^\k.mz).q'^''^^\rnz)) (4.86) 

fc =0 mj =0 

where 


r,(i/-t-l) _ p(i') ^ _ p - ,773 


^(■^+ 1 ) = - m 3 

”^ 3 ) = (-!)' 


q^-) \ / p^'') - 1 + d - /r - nj3 \ A% - q^'''^ ) 

7770 / V d — k — m 3 / 




k W d-k 
m 2 ) \ d — ‘m2 

i — I + 1712 — k 

t - 1 


(4.87) 


vhere the difference between the pole and zero multiplicities in A'M satisfies 


pi'') _ qi'') > 2 



The set of spectral functions that result in the convolution (4.86) differ only in 

the multiplicities and which are determined by the limits of summation 

in (4.86). The ranges of and or equivalently the spectral functions 
are obtained from the order d of the Laguerre cross-spectral function s^iz) and the 
pole and zero multiplicities of X^''^ as follows 

p(u+i) = -f 1,... + d+ 

= 0,1,2, (4.88) 

Now let us consider in detail the entire convolution (4.84). This can be 
obtained by following a set of successive steps that give the value of this convolution 
in an organized algebraic procedure. The value of this convolution is a sum of spectral 
functions X^‘*^ of the form 


where 




(4.90) 


It is clear that the set of zero and pole locations are determined simply from the 
number of convolutions. The value of the convolution is obtained by determining 
the values of pi, qi and Ci based on the orders of the Laguerre spectral functions 
di.d 2 ,“-, d^+i- We can make the following observations: 


1 . If this convolution were performed successively then as a result of the first con¬ 
volutions, Sdj(r) * 5 ^ 3 ( 2 ) we would have pole and zero multiplicities determined 
by a plot similar to that of Fig. 4.2. Each of the resulting functions X^^Xz) is 
then convolved with 5 ^,( 2 ). The pole zero multiplicities in the resulting spec¬ 
tral function ^^(^^(r) are obtained by extending those in Fig. 4.2 using the 
extension map described by Fig. 4.3 for the Laguerre spectral function order 




Figure 4.5: The zero and pole multiplicities in 

m 3 . This results in a range of the pole and zero multiplicities described by a 
new plot similar to that of Fig. 4.4. This procedure is continued until the last 
extension map corresponding to the Laguerre spectral function order m^+i. We 
can obtain the same result directly however by taking the region of pi and qi 
in (4.89) to be as follows: 

• The maximum value of zero multiplicity is dj - 1. 

• The minimum value of zero multiplicity is max[0.di — dj — • • • — d„+i]. 

• The maximum value of pole multiplicity is dj 4- d 2 H- ■ • • + d^+i + 1- 

• the minimum value of pole multiplicity is di + 1 . 

We can also use a plot similar to that in Fig. 4.5. 

2. The coefficients in (4.89) can also be obtained by following a recursive 
operation. Let us assume that after the first v convolutions a set of coefficients 
were computed and associated with the corresponding pi and 9 ,-. For the 
u + 1 *‘ convolution, involving 5^^+,, the extension map in Fig. 4.3 for d „+2 
is used to determine the multiplicity pairs that result from the 

pairs [pW. ^(*')]. Then if from the extension map the v + r‘ convolution results 
in generating from then the coefficient is modified 







to using the multiplier defined in (4.87). This multiplier is now 

written as by to indicate that it corresponds to the multiplicity 

pair generated from the pair after v \ convolutions. Thus the coefficient 
jg updated according to 

^(-+1) = (4.91) 

Finally, for each value of the vectors of indices ot, the Q-spectral function 
can be obtained as follows: 

1 . The multiple convolution of the Laguerre spectral functions which are the en¬ 

tries of one block S" of the matrix S given in (2.106) is obtained as a weighted 
sum of the spectral functions pi,qi) in the form of (4.89). The poles of 

these spectral functions are located at while the zeros are located at 

where i/ equals the sum of the entries of the block M“ of M that corresponds 
to the block S" of 5. 

2 . The range of the pole multiplicities pi, the zero multiplicity qi, and the coef¬ 
ficients Ci in (4.89) are determined as described above from the orders of the 
Laguerre cross-spectral functions in S" for which the corresponding entries in 
M" have non-zero values. 

d. Convolution of 3 q”'^(z;p) with sj^'^^\z:p) 

To obtain the value of the multiple convolution of the Laguerre cross- 
spectral functions in the block S" we first convolve the Laguerre cross-spectral func¬ 
tion 5o(.:; p) a number of times ?7io corresponding to the sum of the entries on the 
principal diagonal of the block M“ of the multiplicity matrix M. The result of this 
self convolution is given by (4.60). We also perform the convolution on the rest 
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of the Laguerre cross-spectral functions (residing on the non-principal diagonals of 
S" the value of this convolution is given by (4.89). Thus the value of the required 
convolution is then given by 

= '£cts'r\z-,p)*X^'Hz;Pi,^i) ( 4 - 92 ) 

t 

The terms under summation in this expression are obtained by convolving ‘’^( 2 ; p) 
with each of the spectral functions X^'^\z-pi,qi). This convolution has the form 

p) . A’W(.i ft, = 2 ^ £_ ?i) 

" " />"*» 2-Kj fc, {v - p^°)iv - p-^'>){v - 

(4.93) 


Since the integrand has only one pole inside the unit circle at location p the value 


of the convolution in (4.93) is easily obtained as 


s<r^\z:p)*X^'^Hz;pi,qi) 


^i(mo + «+l) pTno-^-K-l^qi 

pqi{mo-^K-l) — ^mo+«+l^P» 


(4.94) 


Notice that this convolution results in one term only. Moreover the pole and zero 
multiplicities of do not change and the value of this convolution can be easily 
obtained by just moving the pole and the zero from and to p"*®+'‘+^ and 
^mo+K-i respectiv'ely. More precisely from (4.94) and the definition (4.85) we have 

p) * X^'^Hz-, Pi, qi) = Pi-<ii) (4-95) 
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V. MODELING A NON-GAUSSIAN RANDOM 
PROCESS USING THE WIENER MODEL 


In Chapter IV we develop the procedure for computing the cumulant and polyspec¬ 
trum (of any order) of the output of the Wiener model. These higher-order statistics 
are shown to be functions of the model parameters defined in Chapter III. The first 
such model parameter is the variance of the zero-mean white Gaussian process 
w(n) that drives the system. This parameter specifies the set of Hermite polynomials 
that describe the memoryless nonlinearities represented by section 2 of the model. 
The second parameter, p, defines the pole and zero locations of the transfer functions 
of the linear filters comprising section 1 of the model. The impulse responses of these 
filters are the discrete Laguerre functions and the outputs are zero-mean Gaussian 
processes each with variance cr^. Thus, the auto- and cross-correlation functions of 
the outputs of these filters (and in fact all of the higher-order statistics) are com¬ 
pletely specified by and p. The remaining parameters are the coefficients cct in 
the Q-polynomial representation of the output of the model (3..34). According to our 
definition of the indices represented by the vector a, these parameters are uniquely 
determined once the model dimensions Nl and Nn are specified. Thus for a given set 
of model dimensions the number of the model parameters is finite and the expressions 
of the model output statistics can be obtained. 

In this chapter we demonstrate the application of the Wiener model in repre¬ 
senting some non-Gaussian discrete random processes that fail to meet the conditions 
of linearity. The parameters of the discrete Wiener model are determined such that 
the model output statistics match the modeled process statistics up to some given 
order over a specific region of support. 



A. NONLINEARITY DETECTION IN RANDOM 
PROCESSES 

The statistics of a linear random process either in the time domain (cumulants) 
or in the frequency domain (polyspectra) have characteristics that nonlinear processes 
do not share. A significant body of work has been reported to detect and characterize 
the nonlinearities in time series [4, 33, 5]. 

In [34] Lawrance and Lewis (1987) present a class of stationary random processes 
x(n) for which the autocovariance function, autoregression structure, 

i.e., this function satisfies the p*'* order Yule-Walker autoregression 

- 1) + - 2) + .. • + - p) (5.1) 

where a,- are appropriate real-valued constants. This class of processes includes 
(among others) two subclasses of stationary random processses with mean p that 
have an autoregerssive representation. The first is the class of processes for which 
the linear conditional expectation has the form 

£[x(n) - p\x{7i - l),a-(7? - 2), • • • ,.r(7? - p)] = 

ai(x{n - 1) - //) + a 2 {x{n - 2) - p) + ••• + apix{n - p) - p) (5.2) 

Processes in this class may be nonlinear. The second class, which is a subset of the 
class satisfying (5.2), is that of linear stationary random processes with mean p which 
satisfy 

x{n) - p = ai{x{n -1)- p) + a 2 (x{n - 2) - /i) +- \-ap{x{n-p) - p) + w{n) (5.3) 

where w{n) are independent and identically-distributed (i.i.d.) and the coefficients 
Gi are real-valued constants satisfying the condition that the linear autoregressive 
model represented by (5.3) is stable. 



The order p linear autoregressive residual c(’’)(n) is given by the process x(n) 
after subtracting its best linear least-square predictor in terms of x{n — 1), 
x(n -2),---,x(n -p), i.e., 

c(*')(n) = x(n) —p —ai(x(n —1) —p) —a 2 (x(n —2) —p)- ap{x{n — p) — p) (5.4) 

Lawrance and Lewis [34] based the principle of detecting the nonlinearity in processes 
with autoregressive correlation structure (5.1) on the following: 

1. For processes satisfying (5.1) the linear autoregressive residuals, are un¬ 
correlated. 

2. For linear processes satisfying (5.3) the linear AR residuals are not only 
uncorrelated but they are also indtpcndcnt. 

According to this distinction, the higher-order dependency of the residuals is inves¬ 
tigated by considering the two cross-correlation quantities : 

Corr[(x(n) — //)^. + /)] 

Corr[(x(n) — p)-+ O) ] 

for I = 0, ±1.±2. •••. If a process satisfying (5.1) is linear (i.e satisfying (5.3)) then 
the two cross-correlation functions in (5.5) are zero for all I ^ 0. This test is a useful 
one, but it is based on a particular class of processes defined by (5.1). 

A more general test for nonlinearity was formulated by Rao and Gabr (1980) 
[.35] and later extended by Hinich (1982) [36]. This test is based on the fact that if a 
process x(n) is linear (in mean square), then there exists a sequence h{n) such that 
x{n) is represented by 

OO 

x{n) = ^2 ~ (5.6) 

k=0 
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where «;(n) are i.i.d. and the sequence h{n) is square summable . In this case, given 
the process bispectrum and power spectrum it is easily shown 

that the ratio 

__ (5j) 

5i'^(u>i)5?^(a;2)5i"^(a;i+a;2) 

is a constant independent of u>i and ufj. In practical applications estimates of the 
process bispectrum and power spectrum functions are used to test for its linearity. 

The linearity test of Hinich actually follows from the fact that the k -order 
polyspectrum of a process represented by (5.6) is has the form 

-W2-Wfc-l) (5.8) 

where j3k is a nonzero constant and is the Fourier transform of h{n). In other 

words, linear processes have the characteristic that their polyspectra are factorable: 
this is a generalization of the correlation-based innovations representation of random 
processes. The representation (2.4-3) is possible if the process power spectral density 
function is factorable, i.e., if it satisfies the Paley-Wiener condition. For a rational 
power spectral density, the procedure developed by Wiener [1] can be used to deter¬ 
mine the poles and zeros of the linear system transfer function. The factorization of 
polyspectra however is not as simple. 

Tekalp and Erdem (1989) [13] established necessary and sufficient conditions for 
the existence of a stable linear time-invariant system driven b\ a non-Gaussian k 
order stationary white process such that the system output /c‘^-order polyspectrum 
matches a given one. The procedure is outlined briefly below. 

The A-‘'‘-order complex cepstrum of a k^^-order stationary random process x{n) is 
defined as the multidimensional inverse ’-transform of the logarithm of the complex 
polyspectrum .^2? • ■ • • ^k-i) 

h----- k-i) = {log[5W(zi. .’2. • • •, -'fc-i)]} (5.9) 
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The region of convergence is assumed to contain the fc-dimensional unit circle 
ki| = kal = = kfcl = 1? the integration is performed over that contour. 

For a discrete random process that has the representation (5.6) its A;‘^-order complex 
cepstrum is given by 


fc-i *-i 


hi"' 1 h-i) = f^wHhi hi"' 1 h-i) + 9{~h) H ^») + ^ 9{h) H 


»=1 i=J 

J¥i 


for any i 


(5.10) 


where 6(/) is the unit sample function and 


gil) = Z-^{\og[H(z)]} 


From (2.43) the complex cepstrum of x{n) is also given by 


fc-i fc-i 


ci^Hh-h-'" Jk-i) = logL^fc] Ylih) + 9{-h) n ~ ^‘) + FI (5.12) 


1=1 j=i 


for any i 


Equation (5.12) implies that the region of support of the complex cepstrum of x(n] 
is the union of the regions 72-; defined as 


Tli = {{hihi"'ih-i) ■■ -oo < /i < oo./j = 0 for ; = 1.2. • • •, A' - 1 and j ^ i, 


? = 1,2,---,A- 1 


Hk = {{hi hi"'i h-i) • -oo < /i = /j = • • • = Ik-i < oo} 


(.5.13) 

(5.14) 


Tekalp and Erdem proved that a necessary and sufficient condition for the process 
x(n) to have the linear representation (5.6) is that its complex cepstrum is zero 


outside the region of support 
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If this condition is not satisfied the polyspectrum is not factorable and there does not 
exist a linear time-invariant system which when driven by a white process gives the 
required output polyspectrum. In this case a nonlinear model needs to be considered. 

The likelihood of encountering a factorable polyspectrum has also been studied 
by Tekalp and Erdem [12]. They proved that the subset of factorable polyspectra has 
measure zero in the set of polynomial polyspectra that correspond to finite-support 
cumulants. Thus from a practical point of view the linear non-Gaussian model can 
be used only as an approximation when this is justified. 

B. THE DISCRETE STOCHASTIC PROCESSES 
REPRESENTABLE BY THE WIENER MODEL 

The structure of the discrete Wiener model (in Chapter III) together with the 
analysis of the model’s cumulants and polyspectra (Chapter IV) provide the following 
characterization: 

1. The linear stage of the model, section 1. is represented by a bank of causal, 
stable linear filters whose impulse responses are the discrete Laguerre functions. 
The transfer functions of these filters. Ai(^), are analytical functions with no 
zeros on the unit circle. The outputs of this stage are zero-mean Gaussian 
random variables and their cross-correlation functions are given by the Laguerre 
cross-correlation sequences ri{l). Thus, the complex cross-spectral functions of 
these outputs are given by the Laguerre cross-spectral functions s<i( 0 ) which 
are also analytical functions with no zeros on the unit circle. 

2. The nonlinearities in the model are represented by the Q-polynomials which 
are a wieghted sum of multinomials of the outputs of section 1. Since the model 
output is a sum of products of correlated Gaussian variables, the output cu- 
mulant is computable as a specific sum of products of the Laguerre correlation 
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sequences rj{l). The output polyspectrum is given by a sum of multiple convolu¬ 
tions of the Laguerre complex cross-spectral functions 5 ^( 2 ). This convolution 
was shown to be representable as a sum of spectral functions 
which are also analytical on the multidimensional unit circle. 

While sufficient conditions that a random process must satisfy to be representable 
by the Wiener model are not known, certain necessary conditions are evident from 
the above characterization: 

1. The cumulants of any order of a processes to be representable by the Wiener 
model must asymptotically approach zero as any of the time lag arguments 
increases. This follows because output cumulants of any order of the model are 
obtained as products of the Laguerre cross-correlation sequences ri{l) and all 
these sequences asymptotically approach zero as their lag arguments increase. 

2. bnless the polyspectrum is zero everywhere, as in the case of a Gaussian pro¬ 
cess. a process to be represented by the Wiener model must not be bandlimited. 
a bandlimited process is here defined as one for which at least one polyspec¬ 
trum of order k is zero over a nonzero-volume in the Ar-dimensional subspace 
that defines its region of support. 

3 . The process to be represented by the Wiener model must not be periodic or 
almost periodic [37]. Since the model polyspectrum is represented with spec¬ 
tral functions that have the same form and all have multiple poles at the same 
location, these spectral functions do not have any singularities on the mul¬ 
tidimensional unit circle. Therefore these functions cannot produce any line 
spectral components. 
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C. DISCRETE PROCESS MODELING BY CUMULANT 
MATCHING 

For a given model dimensions Nj, and Ns the output of the Wiener model is 
given by 

= (5.15) 

0!t 

and the k^^-ovdev output cumulant is given by 

Ci^Hl) = E^Ca{l) (5-16) 

a 

where each a vector is a concatenation of the k oii vectors and the summation 
is only over those combinations for which the are non-zero. Let us define the 
vector c as the vector of N coefficients appearing in (5.15) and ordered with indices 
in ascending order, where order for the indices is defined by (4.3). Thus 

c = [ca,-ca3 --*-.cajv]^ 

where here the variables a^.a,. • • • and as represent specific vector of indices of the 
N terms that appear in (5.15). Further define the vector of cumulants 

= lCSo(') ■■■ (5.18) 

where the variables represent the paritcular larger vectors of 

indices a that appear in (5.16). Then (5.16) can be written using matrix notation 
as 

CW{1) = (cW(/))M‘ (5.19) 

where 

c®* = c (8) c 0 • • • 0 c (5.20) 

-- - -^ 

k times 

and 0 denotes the direct (Kronecker) product of vectors. 
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If a set of specific lag values are now chosen, then a vector of these 

values can be formed as 

X^=CWc®* (5.21) 

cWT(ii) 
cWT(/2) 

cWr(/M) 

The entries in each row of the matrix are the Q-polynomial cross-cumulants 
which are functions of the parameter p and the input variance a\. The rows differ in 
their time lag vector. 

Let us observe that although the quantities and p are parameters of the 
system, the theory provides a model using any finite positive real value for cr* and 
any p satisfying \p\ < 1. Although the variance al is essential for ensuring the 
orthogonality of the Hermite polynomials, once its value is chosen, the scale of the 
output can be adjusted by appropriate scaling of the cq;^ parameters. Therefore in 
our modeling, we fix the value of al at 1 and do not attempt to adjust it. 

Although in theory the model could be applied with any fixed value of p with 
magnitude less than 1 , in practice, because the model dimension Nl is finite, the 
value of p has a significant effect on the error involved in matching a given set of 
higher-order statistics. Therefore, when developing the Wiener model for a random 
process, we include both the parameter p and the coefficients coj in spite of the fact 
that there could be some redundancy. 

1. Model Parameter Identification Using the Extended 
Kalman Filter Algorithm 

Given an appropriate stochastic process represented by a set of measurement 
data, we seek to model it as the output of a Wiener discrete nonlinear system. In this 



(5.22) 



case we estimate the process cumulants over defined regions of support and arrange 
these estimates in a vector Xp. The vector Xm in (5.21) is then defined over the 
same region of support. Let us define the modeling error E as 

E = Xp-X„ (5.23) 

We seek to determine the values of the model parameters that minimize the squared 
magnitude E*^E. Since the components of the vector X^ are nonlinear functions of 
the model parameters, a technique for solving a nonlinear optimization problem is 
required. The extended Kalman filter algorithm was used because it was found to 
converge to a solution faster than other techniques tested (e.g., the steepest descend 
method). During the optimization we need to know the cumulant vector X^ and its 
gradient with respect to the model parameters; therefore it is necessary to have an 
efficient method to compute these arrays. In Appendix D we describe a procedure 
that enables one to compute both the cumulant vector and its gradient in a type of 
look-up table that saves signficant computation and storage. The extended Kalman 
filter algorithm, which is used in the optimization, is also described in Appendix D. 

2. Simulation Results 

In this section we present the results of two experimental examples. In 
the first example the procedure developed in this research is applied to model a 
synthetically generated set of data. In the second example the procedure is applied 
to a data set obtained from a biological measurement record. We refer to these 
data sets as the “original data.” In both examples the vector Xp that represents 
the original data cumulant values consists of the estimates made from the data of 
both the second- and third-order cumulants in the nonredundant regions of support. 
The vector X^ representing the model output cumulants is constructed for the same 
regions of support. Note that cumulant values used in this vector are computed using 
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the expressions developed in the dissertation, i.e., they are not estimated from the 
model output. The extended Kalman filter algorithm is then applied to identify the 
model parameters minimizing the squared magnitude of the modeling error. After 
the model parameters are obtained, they are used to generate a data sequence from 
the Wiener model using these parameters, i.e., the model is driven by a computer 
generated white Gaussian sequence to produce data. The cumulants of the generated 
output are then estimated and compared to the cumulants estimated from the original 
data. 

In this section we also compare the generated sequence to the original data 
record. This is repeated for some different model sizes, i.e., some different values of 
Nl and N^. 

a. Synthetically Generated Data 

In this example a discrete \\ iener model was constructed with dimen¬ 
sions Nl = A'n = 2. The poles and zeros of the corresponding transfer functions are 
defined by the parameter p wTich is chosen to be 0.65 in this simulation. Section 3 
in this example is defined by the values of nine coefficients ca not including cq. A 
zero-mean unit-variance Gaussian sequence of 1024 time points was applied to the 
model to generate the random process. The second- and third-order cumulants of 
the random process were estimated as shown in Fig. 5.1. In the following experi¬ 
ments models of different size were used to model the generated data sequence. The 
Laguerre dimension Nl was given values in the range {1.2,3.4,5} while the nonlin¬ 
earity dimension Njf was given the values {2.3.4} for each value of Nl- From the 
results of these experiments we could form the following conclusions: 

1. For all the model dimensions, the estimated value of the parameter p was very 
close to the original value namely in the region [.635. .671]. 
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2. None of the trials resulted in values of the components of the vector c relatively 
close to the original values although in all cases the cumulant matching error 
was acceptably small (< 5%). 

3. Although the cumulant matching v/as acceptable, the plots of the model out¬ 
put sequence showed temporal variations judged to be similar to those of the 
original sequence for the case of quadratic models only. For higher degrees of 
nonlinearity the sequence had temporal variations that appeared to be different 
from the original data. 

The results for {Nl = I, Ns = 2} and {Nl = 2, Ns = 2} are shown here to 
demonstrate the efficiency of the Wiener model in modeling nonlinear data sequences. 
The resulting vectors of cumulants are compared with that of the original data in Fig. 
5.2. The output sequences of the different size models is compared to the original 
data sequence in Fig. 5.3. 

b. Biological Measurement Data 

In this example we attempted to model a data set obtained from a 
record of biological measurements. This data set is one that has been made available 
to researchers in the field of time series anah'sis by the Santa Fe Institute [38]. The 
complete data record contains 34.000 points of the heart beat rate, chest volume, 
and blood oxygen concentration together with the EEG state of a sleeping patient in 
the sleep laboratory of the Beth Israel Hospital in Boston. One segment of the blood 
oxygen concentration data was chosen while the patient was in the same EEG state 
for a sufficiently long period of time (around 2500 points). A plot of the selected 
data segment is shown in Fig. 5.4, while estimates of its second- and third-order 
cumulants are shown in Fig. 5.5. This data set was modeled with the same range 
of dimensions as in the first example. The model parameters were determined such 
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that samples of the second- and third-order cumulants of both the model output and 
the data set match in the least square sense. The results of these trials showed the 
following: 

1. The estimated values of the parameter p were very close together in all the 
trials and approximately equal to 0.95. This indicates that the linear part of 
the model has a low-peiss frequency response. 

2. The vector of coefficients c was considerably diverse among all the trials al¬ 
though the cumulant matching was acceptable. 

3. The resulting data plots were comparable to the original data plot for quadratic 
models but unlike the original data for higher degrees of nonlinearity. 

The cumulant vectors ol the data set and the model output are compared in Fig. 5.6. 
The details of the data and the model outputs for different model sizes are shown in 
Fig. 5.7. 


c. Analysis of Results 

The results of the above examples indicate some characteristics of the 
procedure used for discrete random process modeling. These characteristics can be 
summarized as follows: 

1. The linear part of the model adjusts rapidly to the frequency content of the 
modeled process. The filter bandwidth is properly adapted to the input process 
bandwidth as shown by the estimated value of the parameter p. 

2. The solution obtained for the coefficients of section 3 of the model is not unique. 
This can be explained by the sparseness of the matrices constructed from the 
theoretical output cumulants. The resulting nonlinear equations to be solved for 
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these coefficients can be partitioned into a group of uncoupled sets of nonlinear 
equations, each with a separate subset of these coefficients. Each of these sets 
of equations has more than one solution. One solution from each set provides 
a problem solution that exactly gives the same result as any other combination 
of the different solutions of the different sets. 

3. Over estimation of the model nonlinearity order may give an acceptable cumu- 
lant matching result while the time variations of the data sets may appear to 
be far different. This seems to be illustrative of the fact that for a highly non¬ 
linear model second- and third-order cumulants are insufficient to characterize 
the process. 

The method employed here to match cumulants of a given time series are not the 
only possible methods and the results summarized above are as much dependent on 
the method as they are dependent on the model. We hope that other experiments 
and methods will later be investigated. 
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Figure 5.3: Example 1 data (a): Original data sequence (b): Model output for 
= 1, Nff = 2 (c): Model output for Ni = I, A'n = 2. 



n 


Figure 5.4: Example 2 data : a segment of the data set representing a patient’s blood 
oxygen concentration. 

















Figure 5.5: Example 2 data (a): Estimate of data third-order cumulant (b): Estimate 
of the covariance function. 




Figure 5.6: Cumulant vectors comparison for example 2. “o" = original set 
‘'■x’'=model output;(a): Nl — I.Nn = 2. (b): Nl = 2,Npf = 2. (c): Nr = 














VII. CONCLUSION 


The canonical representation of nonlinear systems developed by Norbert Wiener 
to circumvent the problems of their series-like representation has been used for many 
years in parameter identification and characterization for such nonlinear systems. 
The applications that have been reported in this regard have been primarily for 
continuous time systems and are based on second-order statistics. The kernels rep¬ 
resenting the system are obtained by computing the cross-correlation between the 
model output and its driving input. 

The research presented in this dissertation is focused on the use of the Wiener 
nonlinear system model to represent a general class of discrete stochastic processes 
and, in particular on: 

1. Development of the model for discrete nonlinear systems and discrete random 
processes. 

2. Analysis of the higher-order statistics of the model output, which can represent 
a broad class of discrete random processes. 

3. Application of the theory and the results to model empirical data. 

The Wiener model is comprised of three cascaded stages. The first stage is a 
bank of orthogonal linear filters driven by the the same white Gaussian process. The 
second stage consists of sets of identical memoryless nonlinear modules which form 
orthogonal multinomials of the outputs of the linear stage. The model output is then 
represented as a weighted sum of such multinomials. This output is characterized by 
the summation weights which are the constant parameters of the model final stage. 
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The higher-order statistics are generally very demanding in computation and in 
data size especially if they are involved in nonlinear system applications. The Wiener 
model is shown in this dissertation to be tractible in such applications because of its 
inherent orthogonality properties. Since the model input is white Gaussian, the or¬ 
thogonality of its components results in expressing the auto- and cross-correlation 
functions of the outputs of the linear stage and generally all the model output cumu- 
lants in terms of the cross-correlation sequences between the different order Laguerre 
functions. These cross-correlation sequences depend only on the difference between 
the orders of the correlated Laguerre functions. As a result of this property, a very 
small number of these sequences is needed to formulate the model output statistics 
of any order. Another important characteristic of these functions is that (except for 
the autocorrelation sequence) they are single-sided, which makes most of the terms 
in the general expression of the output statistics identically zero. 

We have exploited the model characteristics in developing a procedure that 
minimizes the computation and storage required to obtain the value of the higher- 
order cumulants of the model output. We have also developed the computation of 
these cumulants such that most of the effort is done off-line only once and stored 
in look-up tables. These tables can be used in any application of this model that 
requires computing the higher order cumulants of its output. 

In the frequency domain we introduced a method that reduces the cost of com¬ 
puting the model output polyspectra. The computation of such multidimensional 
spectral functions involves a large number of integral convolutions. We first proved 
that, similar to the cumulant functions, the size of this computation is greatly reduced 
due to the model structure. Moreover, we showed that such convolution computations 
can be performed algebraically in a structured recursive procedure. 
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In the last chapter of this dissertation we discuss the representability of random 
processes using the Wiener model. The polyspectra of such processes are necessarily 
analytic and not bandlimited. Two examples of modeling processes satisfying these 
conditions are presented. These examples demonstrate the viability of our approach. 
We hope that further experiments of this type, possibly using other methods of 
matching will be conducted in the future. 

The development of this structured model still has areas to be covered and ques¬ 
tions to be answered. In analyzing the characteristic of discrete random processes 
that can be represented using the Wiener model we provided some necessary con¬ 
ditions. Sufficient conditions are yet to be found and the class of those processes 
representable by the Wiener model needs to be defined. 

The development of a procedure to compute the model output polyspectra can 
be done following the results of the time domain development. The algebraic formu¬ 
lation of computing the multiple convolution of Laguerre cross-spectral functions can 
be directlv applied to obtain the power spectral density and higher-order polyspectra 
of the model output. The use of these results to predict model behavior and for 
svnthesis of processes with given polyspectral characteristics has yet to be explored 
however. 

The spectral functions denoted by have multiple poles and zeros 

at locations and respectively. These poles and zeros move towards the 

origin as o increases. This means that as the number of convolutions required to 
represent the polyspectra increases, these spectral functions approach the form 2 “*. 
Since the number of convolutions is proportional to the degree of nonlinearity, in 
highly nonlinear terms the spectral functions can be approximated by 2 "* which 
greatly simplifies the representation. In effect, these terms can be approximated 
by a weighted sum of functions of magnitude one with linear phase. Even if the 
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nonlinearity is not high we can extend the number of Laguerre filters to compensate 
for a reduction in the magnitude of p to justify this approximation. These alternatives 
need to be analyzed and the effect of each on the the model output representation 

needs to be evaluated. 

In summary, the work reported in this dissertation may be just the introduction 
to a whole new body of work yet to be explored. We have demonstrated the feasibilty 
of the Wiener model as a generic representation for a large class of random processes. 
This model expands in a natural way the much more restricted class of linear Gaussian 
processes that are currently so well understood. The utility and practicality of the 
model result from important structure foreseen by Wiener in its development and 
the recent new interest in higher-order statistics that has undergone development in 
only the last decade. Although our initial exploration of this new area leaves man> 
questions unanswered, we hope that the problems posed here will continue to be 
addressed and lead to a significant improvement in our ability to understand and 
model the signals and other data that are encountered in the real world. 
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APPENDIX A 


EXPECTATION OF THE PRODUCT OF 
JOINTLY GAUSSIAN RANDOM VARIABLES 


Let y = [yi,y2,’ ’ ' iVk] be a vector of real zero-mean jointly Gaussian random 
variables. The expectation of the product of these variables raised to some powers 
has the form 

fWs?’•••!'”} (A-1) 

k 

The value of this moment is zero if the sum of the exponents ^ j/,- is odd. Let us 

»=i 

outline a method here for computing this moment when the sum of the exponents is 
even. The value of the moment can be computed by generalizing equation (2.34) to 

= £{ yiyi • • • yiy2y2 • • • yz/ • - ykyk • ■•yk) 

1/1 i/k 

w’here the exponent qj^j^ means that j/jj and yj^ are paired qj^j^ times and the sum¬ 
mation is over all the possible permutations. Since more than one permutation can 
result in the same value of the term under summation (A.2) can be put in the form 

(A.3) 

p 

where in this case the summation is over all the distinct permutations and the coef¬ 
ficient Cp is the number of non-distinct permutations that have the same value. The 
multiplication is over all the distinct pairing of the pj's with each pair raised to a 
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power that equals the number of its repetitions in the pairing configuration. In all 
cases the relation 

ixth *=1 

must hold for each term under the summation. The procedure to find the value of 

this moment begins by constructing the k y. k correlation matrix C 

■(7(1,1) (7(1,2) C(l,3) ••• C(l,fc)‘ 

C(2,l) (7(2,2) C(2,3) ••• (7(2,/:) 

C= C(3,l) (7(3,2) C(3,3) ••• (7(3,/:) (A.5) 

Cik,l) C{k,2) C{k,3) ••• C{k,k) _ 

where 

C{iJ) = SiyiVj} (A.6) 

and the same size matrix M of non-negative integer entries such that 



2?a 

?2 

*3 

' • ik 


ik+i 

‘2ik+2 

ik+3 


M = 

hk+i 

hk+i 

2i2fc+3 

' • izk 


_ i{k-l)k+l 

i(k-i)k+2 

i(,k-l)k+3 * ■ 

• • 2^2 


Note that for clarity in the following discussion it is prefered to use k^ linearly indexed 
variables ij to represent the entries of M rather than variables with dual (row,column) 
indices. W'e also force the diagonal entries to assume even values for reasons that are 
explained in a moment. 

The entries of the multiplicity matrix M are used to determine the exponents 

in (A.3). Since, in general, some of the exponents of the random variables yi 
in (A.2) are greater than one, there are some permutations that result in pairing 
a random varaible t/j with yj more than once. This results in having the quantity 
C{i,j) appearing in the product in (A.3) a number of times equal to the number of 
the resulting identical pairs of ?/,• and y,. The value of the off-diagonal entry M(?, j) 
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is the number of yi paired with the same number of to give The 

diagonal entry M(j, j) in row and column j is the number of the elements taken from 
y^ that are paired together and the result gives C{j,j) * . Therefore the values of 
the diagonal entries must be even. 

At this point it is required to compute the number of permutations to perform 
the pairing for a specific configuration of the entries of the matrix M. We start from 
the first row and examine all the possible values that the entries of this row can 
assume. Then we examine the other rows moving downwards in order. When a row 
is examined its entries are examined in order from left to right. Therefore starting 
with the upper left entry M(l,l) = 2?! for ii = 0,1,2, • • •, [yj the number of 
permutations is obtained by computing the number of ways of choosing 2?! elements 
out of 1^1 multiplied by the number of pairing permutations of 2i\ elements. This 
yields 


per 


2ti!(i/i-2?a)! ii!2*i 

- 2?i)!2“ 


(A.8) 


Since the sum of the elements in the first row of M equals then for each value of 
7 i the value of M(l,2) = < Vi - 2?'i. Also the sum of the elements in the second 

column must equal 1/2 which means that 22 < 1 ^ 2 - Therefore the multiplicity matrix 
entry M(l,2) equals ^ such that ?2 = 0.1.2,• • •,min{i/i - 2i-^,V2}. The number of 
permutations in this case is given by the number of ways of choosing elements out 
of - 2/1 multiplied by the number of taking the same number /2 of the y 2 out of 1/2 
then multiplied by the value corresponding to the number of pairing permutations. 
That is, 






P2- 


/2!(i/i - 2/1 -/2)! 


X l2\ 


1.51 



_ _ (t/1 - 2ii)\u2^. _ 

W-{vi - 2ii — - * 2 )! 

Similarly the other entries of the first row of M follow as 
Arfij) “ £j[=i M(l, k))\ _ i 


(A.9) 




- z^k=i •'3- ^ .1 


and the total number of permutations of a specific configuration of the elements of 
the first row equals the product of the numbers of the above permutations,i-c., 


i/i!i/ 2! • • • i^fc! 


„(1)__ (A.ll] 

~ 2^Hy\i2\ • • • 4!(l^l - 221 - 22-jfc)!('^2 “ * 2 )! ■ ' ’ {Vk - 4)! 

For the second row we start with M(2,1) = ik+i 

ifc+i = 0,1,2, • • • ,min{( 2 /i - 2ii - -4), ( 1^2 - * 2 )} 


The number of possible permutations is 
(//i-2?i-22- 


( 2^2 - 22 )'- 


P *’’ 2*+ i !(2^1 - 2 ?i - l2 -2 fe -2 fc + l )' 2 fc + l !(2^2 - 22 - 2 fc + l )! 


_ _ (i/i - 2?i - ?2-?fc)!(2^2 - 22)1 _ 12) 

2fc+l!(2^1 “ 2ii — 22 • ’ • — ik — 2fe+l)!(2^2 - 22 — 2fc+i)! 

For the entry M(2,2) = 2 ?fc +2 with 2 k +2 = 0.1,2. • • ♦, J the number of 


permutations is 


(U2 — 22 ~ 2fc+i)! 


271^,! 


. r ( 22 ) _ _ {‘'2 -•'2 ^ 

' “ 27 fc + 2 !( 2^2 - 22 - 2 fc+l - 24 + 2 )! 4 + 2 ! 2 ‘*+> 

_ _ (t/2 — 72 ~ 2fc+l)- _ 

2’*+*7fc+2K272 — 72 ~ 2fc+l — 27fc+2)l 


(A.13) 


Proceeding with the rest of the entries of the matrix M. we can examine all the 
possible values of the entries and compute the number of permutations of each con¬ 
figuration noting that in all cases the sum of the entries of both the row and the 
column must equal Uj. The coefficient Cp in (A..3) is given by 


cp—n n 

t=i i=i 


(A.h; 
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Due to the cancellation of of terms that appear in both the denominator and the 
nominator, the value of the moment is given by 

31 n .n (MU., j,))! 


/ A 1 


where the summation is over all possible configurations of the matrix M that satisfy 
the above mentioned condition on the sum of the entries along the rows and the 
columns. We can write this condition as 

3 

E M(z-,;) = !/,• (A.16) 

t 

The first product corresponds to the product of the diagonal entries of the correlation 
matrix C (the autocorrelation) while the second corresponds to the off-diagonal ones 
(the cross-correlation). 
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APPENDIX B 


THE WIENER G-FUNCTIONALS FOR 
DISCRETE NONLINEAR SYSTEMS 

A. THE RELATION BETWEEN THE DERIVED WIENER 
KERNELS AND THE LEADING KERNEL 

Since the Wiener Functional ^p[^p,5p-i(p), • * •, <7o(p); is required to be or¬ 

thogonal to any Volterra kernel of order less than p when the input is Gaussian,(3.14), 
the set of p equations 

£|Wi[/ii;ic(n)]^p[ 5 p.Srp-i(p),---, 5 o(p);«^(«)]} = 0 for z = 0,1,2, • • • ,p - 1 (B.l) 

are used to determine the relation between the derived kernels and the leading kernel. 
Equation (B.l) can be put in the form 

S \T-Ci\hi] Le(n)]^p[pp,Pp_i(p), • • • ,Po(p); ^ Gj[ 5 j(p); w(n)]| 

= ^5{?ii[/?.-;ie(n)]G,-[p,(p);uKn)]} = 0 for i = 0,1,2, • • • ,p - 1 (B.2) 

3=0 

where Gj[pj(p); w{n)] is the homogeneous functional corresponding to the leading 
kernel Pp for j = p or one of the derived kernels otherwise. Each of the terms under 
the summation has the form 

{ OO OO oo I 

* E hiikiA' 2 . • • •, ki)gj{p){ki+i. • • •, ki+j)w{n - fci) • • • w{n - ki+j) ^ 

ki=0k2=0 J 

oo oo oo 

= ^ ‘ • • • • ■ ^*+2-■ ■ ■ • ^>+i) 

fcj =0^2=0 

x£{u>(” ~ ki)u'{n — k 2 )---w(n — ki+j)} (B.3) 
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in which the expectation of the product of Gaussian variables is involved. The result 
is that the value of this term is zero if the number of the Gaussian variables i + j is 
odd. This means that if the Volterra kernel hi has even order i, then (B.2) will have 
terms that correspond to ^-kernels with even orders only. Likewise, if hi is odd (B.2) 
will contain only odd order ^-kernels. Since p may be even or odd, Pp may appear in 
either collection of terms. Therefore, for simplicity in the following development p is 
assumed to be odd but the result also applies to the case of p even. 

Starting with the expectation of the product of the G-functional with the zeroth 

order Volterra kernel we have 

{ oo oo 

Mo(p) + E S hog 2 {p)ifcuk 2 )w{n - ki)w{n - kz) 

*,=Ofcj=0 

+ E S 12 kogi(p){ki, k 2 , ks, k 4 )w{n - ki)w{n - k 2 )w{n - k 3 )w{n - k^) 

fc,=0 fcj= 0 jltj= 0 fc 4=0 

+ •••+£••• £ hogp-^p){ki,---,kp-i)w(n-ki)---w{n-kp-i)\=0 

*1=0 fc,,_i=0 J 

(B.4) 


Then applying the properties of the expectation of the product of Gaussian variables 
and the symmetry properties of the Volterra and Wiener kernels with respect to their 


arguments produces 


ho 


go(p) + <3^0 12 92(p)ik, k) + 5^ X] 94 {p)iki,ki, k 2 , kz) + 

k =0 * 1=0 * 2=0 


+ 


(p~ ^)- c,(p-i) Y 

(^)!2(*^) “ *‘^c 


12 9p-i(p){ki,ki, k2, kzr’-i ^i>^> 

*j^=0 


= 0 


(B.5) 


Since this is required to be true for any ho, the zeroth order kernel relation is 
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5 o(p) + <^o S 92 {p){k, k) + 3 <t^ 53 E Hp){ki, ki, ^2, ^2) + • • • 

k =0 * 1=0 * 2=0 

A ■ ■'. *=• • • • ■ =0 

(B.6) 

Now we proceed through the Volterra kernels of even order. Letting t = 2 in (B.2) 
we have 


« E E'‘ 2{ki, k2)go{p)'w{'^ - ki)w{n - Ara) 

[ *1=0 *2=0 
00 00 00 00 

+ E E E E* 2 {ki,k 2 )g 2 (p){k 3 - k4)w{n - ki)w{n - k2)w{n - k3)w{n - k^) 

*1=0 *2=0 *3=0 *4=0 
00 00 00 

+ E E-E‘ 2(^1, k2)g4{p)ik3. A'4, /^s, A: 6 )ic(r? - ki)w{n - ^2) • ■‘w{n - ke) 

fci =0^2=0 fce=0 


00 00 


+ ^ ••• ^"2)5p-l(p)(^"3* ^4* * * • • “ ^’l) ’ ‘ ^p+l) ^ —0 

^1=0 fcy4.1=0 

(B.7) 


Now, let us divide the pairing permutations into two sets the first of which is such 
that u'{n — k'l) and u-(7j — ^- 2 ) are paired together and the second such that they are 
not. We have 
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*=0 


fl'o(p) + 5 ^ 52 (p)(^i? ^i) + ^ 9 *{p)i^^^^ ^2, ^2) 

fcj=0 *i=Ofcj=0 






+ 2<t^ 52 ^ 2 (^ 1 ? ^ 2 ) 

ik,=:0fc2=0 


52(p)(^1, ^2) + 6 52 ^ 4 (p)(^l 5 ^2, ^ 3 , ^3) + • ' • 


fcs=0 


'*’ ^ ot/i-3Mof^ 52 •“ 52 9p-i{p)iki',h,h,h-'-’,ke±y_.ke:^) 

2 ^ =0 


= 0 


(B.8) 


From (B.6) the first term (in brackets) is identically zero. Therefore, the second 
term must equal zero for all values of h)', this provides the second even-kernel 

relation 


92(p){ki, h) + 6 52 9i{p)iki^h, k.k) + 


k=0 


+< 7 ; 


P-3. 


(p- 1)! 


^52"' 52 9p-'i{p)iki'!ki.k3,k3.‘• • ,k^,k^) — 0 


2!(H^)!2(^)*f:'o .4r=o 


(B.9) 


Let us take this one step further by letting ?' = 4 in (B.3) and dividing the 
pairing permutation into three sets. In the first set all of the time arguments in 
h 4 {ki,k 2 ,k 3 ,ki) are paired together and are not paired with any of the p-kernel 
arguments. In the second set only two of the arguments are paired together and each 
one of the other two is paired wdth one argument of the p-kernels, if possible. In the 
third set each of the arguments is paired with one of the arguments of the p-kernels, 
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if possible. By doing this, the value of the expectation we obtain is 


11 ^ 4 (^ 15 ^ 15 ^ 2 ,^ 2 ) 9o{p) +(^l^92{p){k,k) 

ik,=Ofci=0 [ *=° 

OQ 00 

“I" 3(7^ ^ ^ ^ !! 54(p)(^35 ^3, ^4, ^ 4 ) 

fc3=0fc4=0 

+ (p |)- ^ ... ^ ^P_l(p)(fc35^3,^45^45**'5^2±l5^e^) 

(V)!2^ * U,=o *tt,=o 


00 00 00 


+ 12(7® E E i: ^ 3 ) 1 52 (^ 11 ^ 21 ) 

fcl=0 ^2=0 ^3=0 

^ , X f>-3 (p “■ ^)- (p "" 

+ 6(7* ^ 54(p)(fci,^'25fc4,/^4) + '-- + «^o 2!(t:3)!2(^)2!(E^)!2('^) 


13 53 ^p-l(p)(^*l'^2-^'4, ^45 " • 5 ^'t^) 


fc4=0 *E^=0 


00 00 00 00 


+ 24^0 E E E E ^ 4 ( A'l, A'2 • ^"3 ? ^4 ) j P4(p) ( ^’1 • ^"2 ‘ ^3 • ^4 ) ”l~ 

A;x=0 ^2=0 ^3=0 ^4=0 


(p- 1)! 


° 4!(H^)!2(^).^o *4r=o 


Ot) 

13 ■ ■ ‘ 13 S'p-ICp)!^'!-^'2> ^ 3 , ^' 4 > ^' 5 -^’ 5 , • • • 5 ^'E ^5 


(B.IO) 


This is similar to the previous case; the first and the second terms (in brackets) are 
identically zero from (B.6) and (B.9). which implies that the third term must be 
equal to zero for all h^ih^ki^h-h). This leads to the third even-kernel relation 


P4(p)(/ri, ^2, kz, ki) + Ibal ^3 96{p){ki^k2. k^. k^. k^, fcs) + • • • 

*8=0 


+ 


E ••• £ gp-i{p){ki.k2,kz.k4.k,.,k5.->-,k^.,kt^)-0 


4!(^)!2("^)*t^o k^= 


(B.ll) 


By proceeding in the same manner up to i = p - 1 we obtain 


^ 2 ^ terms equal to zero from above equationsj + 

00 oo 

*J=0 *y-l=0 


, kp-i) — 0 


(B.12) 


which implies that 

53 hp_i(A:i,^2,---5 Vi)S'p-i(p)(^i>^2,---A-i) = 0 (B.13) 

*,- 1=0 

for all hp.i{ki, hr", ^p-i)- This relation is satisfied only if 

5p-i(p)(^i,^2,---,Vi) = 0 (B.14) 

which is the last even-kernel relation we need. If we substitute (B.14) in the relation 
just preceding it we obtain 

9p-3{p)ih- h-, • • • 1 ^p-s) = 0 (B.15) 

By continuing this backward substitution in all the above even-kernel relations (B.6)- 
(B.14) we find that all the even-order derived kernels are equal to zero . 

To obtain the relation between the odd-order kernels we perform the expectation 
of the product of the G-functional and odd-order Volterra functionals. Making use of 
the properties of the average of the product of Gaussian variables and the symmetry 
properties of both the Volterra and the Wiener kernels, we can derive the following 
relations. Starting with the first order Volterra kernel we obtain 


oo oo 

<^l Z) * 1 (^ 1 ) 9Hp){h) + 3<t* 53 93{p){h, 

* 1=0 * 2=0 

eo oo 

+ 15a* Z Z 9s(p){h,k2, ^2, h, /Ja) + • • ’ 

kt=0ki=Q 

(B.16) 

Since this must hold for all hi{k) then the first odd-kernel relation is 

OO 0000 

9iip){ki) + 3(7^ 53 93{p)ih, h, ki) + Iba^ Z Z 96(p){h, h, k^, kz, kz) -f • • • 

n >->n —n 


*2=0 *j=0fc,=0 

+P ••• Z 9piki,k2,k2,kz,kz.,-■ ‘ ■,kt:^,k2^) 

(E|i)!2(V) ’ * 


(B.17) 


Using the the third-order Volterra kernel we have 

0000 OO 

3(7^ Z Z kz(ki, ki. k2) 5l(p)(^'2) + 3(7^ 53 5'3(p)(^'2, ^3, ^ 3 ) 

*,= 0 * 2=0 * 3=0 

00 00 

+ 15^0 Z Z 9s(v){k2, kz, k'z. kjj, kj) + • • • 

*,= 0 * 4=0 

+ E ■■■ £ g,(k,.h,h.k,.k„-.k^,k 


/C3=0 fcp43 :=0 


43 ^ At y43 
2 2 


4 3!(7® £1:1:* 3( A’l, A '2 , A- 3 ) 1 ^3(p) (ki.kz-kz) 

* i = 0 * 2 = 0*,=0 

00 

4 60(7^ 53 55(p)(^'1 • ^'23 ^'33 ^’4- ^' 4 ) 4 • • • 

* 4=0 

+ TTT^4‘''£--' E gp{h.h.h.kt.k,.---.,k^,ke^) 

(^)!2^ 2 > k^=o *j^=o " ’ 


(B.18) 
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Since from (B.17) the first term (in brackets) is identically zero, the second term is 
required to be equal to zero for all h^. This implies that 

OO 

9z { p ){ h , h , h ) + 60<7* Y , ^2’ ^3’ ^4) + * * * 

*4=0 

(«f^)!2( » ) fc4=o k^=o 

(B.19) 


This operation is continued until the expectation of the product of the G-functional 
and the last two Volterra kernels hp-t and /ip _2 is performed to yield the last two 
odd-kernel relations: 

9p-4{p)iku k2, • • •, fcp- 4 ) + 


00 00 


+ 


2!(p- 4)!2= 


<^0 X X 9pih.k2.‘--,kp^z,kp.^,kp.2,kp.2) = 0 (B.20) 


fc,_s=0fc,_2=0 


and 


gp-2(p){k\-k2,--- ,kp-2) + (p ^ 2)!• ^’2: • • • ’ Vi ’ ^ (B.21) 


This relation is equivalent to 


-2\ 1 


9p-2{p){k\- ki, - •" ^ kp- 2 ) — 




Now substituting this in (B.20) and rearranging terms we find 


5p-4(p)(^li ^2; • • • 1 kp- 4 ) — 




/ .,2 \ 1 00 00 


2!(p-4)!''“'■ l!(p-2)!V2 


) X X 9piki',k2-.’‘ ’ ■,kp-3,kp-z-,kp-2ikp-2) 

fcy-> 3=0 fcy- 2=0 


OO OO 




OO OO 


7^ 1 ( T S ^ ■ ■ ■ ’ 

2!(p-4)!V2 / fc^_,=ofc,_,=o 


(B.23) 
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Finally substituting this in the preceding relations we obtain 
Sp-6{p)ih, h,"-, kp-e) = 

p\ oo oo eo 

~ 3!(p — 6)! j ^ ^ ^ 9pikiik2',’" ikp-s^kp-i,’• • ,kp-3,kp-3) 


93ip){h,k2,k3) = 
(-1)*T^- ^ - 

^ ^ (^)!(3)! 

9i{p)iki) = 






■ Y. 9p{ki,k2,k3,k4,k4,-‘^,kt^,kt:±) 

f ^ 3 2 

fc^=0 

OO 

■ Y 9piki,k2,k2,---,ke:::i,kt^) (B.24) 

k^-O 


Similar results can be obtained if p is assumed to have an even value. The relation 
between the even-order derived kernels and the leading kernel is identical to those 
obtained for odd p. 


B. AVERAGE OF THE PRODUCT OF TWO DISCRETE 
G-FUNCTIONALS 

In this section we present the proof of the orthogonality of the discrete G- 
functionals, i.e., we show that 

^ {0pA9pi-,'>J^in)]Gp,[fp,;w{n)]} (B.25) 

is zero for pi ^ P 2 - We show how to develop an expression for the value of the average 
for Pi = p 2 = p and give that expression. 

From the definition of the G-functional. Gp[ 9 p'-w{n)] is orthogonal to any ho¬ 
mogeneous functional of degree less than p. When pi ^ p 2 (assuming pi > P 2 ) we 
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expand ^p,[/p,; t«(n)] using (3.14) and substitute its value in the product to obtain 

fl?J 1 

E {Qj^\9pi\w{'n)]Gn[fv2'-i “’(”)]} = ^ j XI [fl'pi; «^(”)]Gpj-2m[/pj-2m; I 

L?J 

= E{Qp^[gpi]w{n)]Gpt-2m[fp3-2m',w{'^)]} 

m=0 

(B.26) 

Each term in the summation is the expectation of the product of the G-functional 
and a homogeneous functional with degree less than p\. Thus from the definition of 
the G-functional, each term in the summation is equal to zero. This proves the first 
part of the property. 

When Pi = P 2 = P all the terms under summation in (B.26) except the first (for 
m = 0) are equal to zero. In this case (B.26) reduces to 

E {gp[gp\w{n)]gp[fp-win)]} = E {gp[gp:w{n)]Gj,[fp:uir,)]} (B.27) 

Now for the sake of clarity let us assume that p is odd and develop the value of the 
expectation. In this case 

^ {Oplgp- ^’(")]Gp[/p;«’(”)]} 

= J ••• S/p(^'i,-",^p)^'(" - ~ 

* 1=0 kf =0 

OO 

fc,»+i=o 

5!^ 53(p)(^>+i* ^V+2^ ” kp^i)w{n — kp^2)^{'^ ^p+s) 

-f •••-!- 53 13 ^p(^’p-i-i’’’'’^2p)^’(^ ~ ^p+i)■ ■ ■ — ^2p) (B.28) 

*,+ 1=0 * 2,=0 _ ^ 

Since the next steps are algebraically tedious but otherwise straightforward let 

us just give an outline of those steps. By switching the order of the expectation and 


164 



summation we involve the expectation of the product of Gaussian variables. Then 
the development is continued in the same way as we derived the relation between the 
kernels. The pairing permutation is divided into sets. In the first set only one time 
argument of the fp kernel is set equal to one argument in each of the g kernels and 
the rest are paired together. In the second set each of three time arguments of the fp 
kernel is set equal to one of the arguments of the g kernels, if possible, and the rest 
are paired together. This continues until in the last set each of the arguments of the 
fp kernel is set equal to one argument of the g kernels. Since for the last set there is 
not enough arguments except in the leading kernel gp we have 
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E ••• E fp{h,k2,k2r".f:^L,kt^) 

\ * 2 " ' /*^ ^ Jkj “0 


X Ui(p)(Ari) + 3cr* 5^ ^3(p)(*i,^«^+n^e±i+i) 


/ _ 1 \| 00 00 

+ -1~ P ^p-i^,gi> ^^o~^ 51 •••5Z5p(^l5^£±l+15^E^+r 

+ /' S • • • £ fpiki,k2, kz, ki, • ^ ■ ,ke±i,ke±t) 

(^)!2 2 fci=o fcj+s 


X i?3(p)(^l,^'2>^3) +60 <To ^ 5S(p)(^'l,^2,^3,^'m+l4'E^+l) 

*^+l 


-f ...-I— _^' ,-T<^o ^ 5 Z • •''^9piki-,k2,kz, ^ + ^yki±s_^_^,-■ • ,kp,kp) 


+ ^ ' * * ' ^P~2^ ^p~l? ^p-l) 

fci=0 fcp^2 

p! * 

^ 9p-2{p){kl, • • • , ^p- 2 ) + _ 9)|9l^o XI 5p(^’l - * ' * ? ^p~2, ^p-li ^p-l) 

+ ■ ■ ■ 53 fp{kii kz, • • •, kp)gp{ki, kz, • * •, /jp) (B.29) 

Jbi=0 fcp 

All the sum of terms between brackets were shown to be equal to zero in (B.6)-(B.22). 
Therefore 

^ {0p[9p- M'^)]0p{fp'^ ^‘^’(”)]} = 53 ’ • • 53 ^ 2 - • • •, kp)gpiki. kz.-”, kp) 

ki =0 kp 

(B.30) 

An identical result is obtained for even values of p. 
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APPENDIX C 


COMPUTATION OF THE NON-ZERO 
VALUED Q-POLYNOMIAL CROSS-MOMENT 


The general expression for Q-polynomial cross-moment is given by (4.23), which 
is repeated here for convenience; 




L%-J L^J 




E[M]=Q:-2m2 3 ji=l 


mi=:0m3=0 mfc=:0 

1 ‘■y (c.o,. (C.U., 




i2=ii+i 




(C.l) 


The goal here is to develop conditions for those terms in the expression which are zero 
or sum to zero since those terms will not need to be calculated. The development of 
sufficient conditions for the general case is very lengthy and requires more space than 
is feasible in this exposition. Instead, let us demonstrate these cases to derive the 
sufficient conditions for the example given in Chapter IV and show’ that the result 
can be generalized. 

We begin by specializing (C.l) to the case oi k = 2 and rearranging the terms 
for our example to write 
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-E 


mi 


^omi!2»"* .^0 .^0 .^0 ii!i2!23!i4!2^‘ ^o ”^4-2”* ^io!2‘» 


t4=aio -2mi -2ti -*a~ts >0 


I ®11 t 

V (-1)"' V “‘v’‘ v (-^r* 1 

7 Tl 2 * 2 ^^ /»_ t /» _ t o ^lOtR 


m2=0 


H:=0 16=0 

t7=aii —2m2—2i5—ie>0 


i6!26!27!2**' ^0 m5!2”‘» ^ll!2‘“ 


B 


^ (-ir^ ry(/) ^ (-ir« 1 

^ „^0 .^0 i8!?9!2‘‘ ^o ”^6!2"** ii2!2‘» 


(C.2) 


t9=ai2 — 2m3’-2ts >0 A 

In this arrangement the values of the rrij and the ij are interdependent due to the 
condition (4.19). Therefore when we rearrange the summation the upper limits des¬ 
ignated by hJ 2 ,i 3 ,U,ts and te must be determined so that none of the ij assumes 
a negative value and results in a rejected configuration of the matrix M. When we 
consider these modified summations in detail there are many well defined cases in 
which terms clearly sum to zero. As we show later, eliminating these cases avoids 
the unnecessary computations which represent most of the terms in the expression. 


Starting with the innermost summation in (C.2) marked A, let us proceed to 
determine the upper limit. From the condition specified by the last row of M in 
(4.19) we have 


*12 


^22 ~ *4 ~ *7 — *9 
2 


— me 


(C.3) 


The following restrictions need to be imposed on the relations between the summation 
parameters: 


1. The quantity 032 — *4 — *7 — *9 must have an even value so that is an integer. 
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2. Since 14,^7 and *9 are non-negative integers and sum to a value less than 022 


Q22 — ^4 ~ ^7 ~ ^8 ^ Q22 


(C.4) 


3. Since ii 2 must be non-negative 


me < 


<^22 “ ^4 ~ ^7 ““ ^9 


(C. 6 ) 


This means that we must set the upper limit of the summation te to the value 
on the right of (C. 5 ) to prevent in from assuming negative values, i.e., 


<>22 — 24 “ *7 ~ *9 


(C. 6 ) 


With these conditions the term A in the expression becomes 


^0 - Tne)! 2‘«<6! 


(C.7) 


which is zero for all non-zero values of te- This means that for some permutations the 
corresponding terms in the cumulant expression are summations of quantities that 
have a common factor of zero. Then to save the wasted effort of computing these 
terms we set the relation between the summation parameters to avoid these cases. 


This means that we let = 0, i.e.. 


022 — 24 — 27 — 29 


(C. 8 ) 


which implies (using C.3) that 


2 x2 = -me 


Since both 212 and me must be non-negative integers the result is 


(C.9) 


me = 2x2 = 0 


(C.IO) 


and the term A in (C.2) equals one. This important relation means that because me 
equals zero we consider only the first term of the Hermite polynomial 7 / 2 ( 2 / 2(22 -f /)). 
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Because ii 2 equals zero, the pairing of the yi does not contain the value 
£{y 2 in + l)y 2 {n + /)} = ro( 0 ) as a multiplier. 

Now, from (C. 8 ) and the relation between ig and zg given in the third row of 


(4.19) we have 


*8 = 


0-12 — 2m3 — ig 


ai2 — 022 + *4 + *7 


— m 3 


(C.ll) 


and the summation over tg disappears because ig is forced to have a specific value. 
In addition, the upper limit of summation over m 3 must be restricted to the cases 
that have non-negative integer values of zg. Therefore, letting 

012 — 022 + *4 + ^7 


the last line in (C. 2 ) starting with the summation over zg is 


“ (-ir^ _ 

^0 m3!2”‘> {h - mzy.{Q22 - U - 77)!2‘*"”*’ 


ra»-i4-»T(/) « (-ipta! 

(022 - U - ? 7 )!< 3 ! 2 ‘^ „“o - ”^3)! 


022 -»4 “-*7 


(/) (1-1)*’ 


(022 — *4 ~ * 7 )! < 3 - 


(C. 12 ) 


Similar to the above discussion this term equals zero for all non-zero values of < 3 . 
Hence the computational effort is avoided if we let tz = 0, which results in 


m 3 = Zo = 0 


(C.13) 


012 = <^22 — Ia ~ 1-7 


(C.14) 


The entire term starting with the summation on m 3 (last line of (C. 2 )) then equals 


rni) 


(C.15) 
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The term marked B in the expression (C.2) is treated in the same way to obtain 


1 X (-Ipts! 

2 ‘»t5? ^0 - ms)! 


(C.16) 


where 


0!2i — is — H 


(C.17) 


This has a non-zero value for ts = 0, resulting in 


ms = zii = 0 


(C.18) 


?6 — ^21 ~ *3 


(C.19) 


which implies that ie has a specific value for each value of ? 3 , and the summation 
over ig disappears. From (C.14) ?7 has a specific value given by 


?7 = Q22 — Oi2 — 24 


(C.20) 


and the condition given in the second row of (4.19) implies that 


h = On — 2m2 — 22^ - ig 


(C.21) 


Therefore from the above 


Oil "I" Oi2 — 021 ~ O22 + 23 ii 

H = ^-m2 


(C.22) 


and the double summation over and ig disappears. The term starting with the 
summation over m 2 (second to last line in (C.2)) reduces to 


,tt21 / 7\„tti2+a22”ttl2~»4 i 


r“--‘>(/)r? 


[-ir^ts! 


(021 — f 3 )!(Ol 2 + O22 — O12 — 24)!t2!2‘’ mj=0 — m2)! 


(C.23) 


where 


On -h O12 — 021 — Q22 + is + i^ 


(C.24) 
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As in the previous cases (C.23) sums to the binomial term (1 1) * which requires 

that <2 = 0 results in the relations 


m2 = is = 0 


(C.25) 


*4 + ia = <^21 <^22 ~ <^11 ~ *^12 

Consequently the term represented by the second line of (C.2) becomes 

(Q21 — i3)!(<^22 — Q'12 + *4)! 

In a similar manner the term marked C in (C.2) is put in the form 

1 X {-ir*u'- 

^0 - ”^4)! 


(C.26) 


(C.27) 


(C.28) 


where 


O 20 — *2 


(C.29) 


and summed to the binomial (1 - 1)‘^ Then it is required that <4 = 0 which yields 


= iin — 0 


77?4 = 7l0 


(C.30) 


Then from (C.26) and the relation 


74 = Oio — 27771 — 27 i — 72 — 73 


(C.3i; 


we obtain 


OlO + Oil + Qi2 — O 20 — 021 — O 22 _ 

7l = -- 


72 = O 20 


?4 = Oil + Q 22 ~ C712 “ <^21 ~ *3 


(C.32) 


To guarantee that ?i is a non-negative integer, we must have 


^ QlO -f Oil -f- O12 — O20 — O21 — O22 ^ ^ 


(C.. 33 ) 
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The total term then becomes 





1. The Q-polynomial cross-moment is computed by computing the expectation of 
the product of only the first terms in the Hermite polynomials in the expansion. 
This corresponds to setting m = 0 in (4.23). 

2. This expectation is taken in a special way such that the pairing permutations 
do not include those that result in one or more autocorrelation functions of any 
of the outputs at zero lag, ro(0). This corresponds to setting the values of the 
diagonal entries of the multiplicity matrix M equal to zero. 
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APPENDIX D 


COMPUTATIONAL METHODS AND 
OPTIMIZATION USING THE EXTENDED 

KALMAN FILTER 


In Chapter V we are concerned with finding a model, through an optimization 
procedure, that matches the a set of cumulant values computed from data. Since 
the optimization is an iterative procedure, it is necessary to have an efficient way 
of to compute the cumulants in (5.21). Also, since the optimization uses gradient 
information, it is necessary to have an efficient wa\’ to compute the derivatives of the 
cumulants with respect to the model parameters. The analysis presented in Chapter 
IV provides the means for this efficient computation. 

A. PROCEDURE FOR COMPUTING THE CUMULANT 
VECTOR 

Since the Q-polynomial cumulants are computed directly from the vectors of 
indices, the structure of the arrays in (5.21) are the same for all the models that have 
the same dimensions AT and Nn- Therefore for every pair of dimensions [Nl,Nn] a 
table-like array of numerical values can be computed and stored in a library. This ar¬ 
ray contains all the information needed in (4.21). Each row in this array corresponds 
to one term in (4.21) and has 3 -|- k{k - l)/2 fields defined as follows: 

1 . The first field has the index i of the k coefficients Ca- i.e., it specifies one 


element in c®*'. 



2. The second field is the power of cr*, which equals half the sum of all the com¬ 
ponents of the vectors of indices a. 

3. The third field is a numerical value that results from dividing the product of 
the factorials of the components of the vectors of indices by the factorials of 
the nonzero entries of the multiplicity matrix M. 

4. Each of the last k{k — l)/2 fields has size Ni -|-1 and corresponds to one block 
R(/) in the matrix Cy(f) in (4.17) with / either one of the components of 1 or 
a difference between two of them. The entries of each of these fields are the 
multiplicity of the corresponding Laguerre cross-correlation sequence rd(/). 

According to (5.21) this array would be A’* long where the number of model coef¬ 
ficients, N, increases rapidly with the model dimensions Nl and N^. This means that 
a slight change in the model dimensions or in the cumulant order would increase the 
computational cost dramatically. Also to compute the value of the cross-correlation 
sequence ri{l) every time it is needed using (3.96) is very costly, especially when 
the parameter p is changed as it is in every iteration of the optimization procedure. 
Fortunately, as we have seen, the model has structure which dramatically reduces 
the computational cost. Recall that; 

1 . The cumulants of the Q-polynomials are identically zero for most of the com¬ 
binations of the vectors of indices a. This can be detected by running a simple 
test on a to determine whether the corresponding term is identically zero or 
not. Using this test, presented in Chapter IV, greatly reduces the length of the 
stored array. This results in a significant saving of both storage and computa¬ 
tion. Tables D.l and D.2 show the actual length of this array compared to the 
quantity for the cases of second and third-order cumulants respectively. 
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1 

2 

3 

4 

4 

13 

29 

54 

T 

25 

81 

196 

9 

45 

145 

370 

9 

81 

361 

1156 

16 

116 

516 

1741 

16 

196 ' 

1156 

4761 

25 

250 



25 

400 



36 

477 



36 





TABLE D.l: The ratio of the non-zero terms in the second-order cumulant function 


Al/An 

1 

2 

3 

4 



42 

225 

824 

1 


125 

729 

2774 

2 


213 

1731 

10060 


729 

6859 

39304 



759 

8746n 

8746 

3 


?r44 

39304 

328509 



2171 



4 


8000 





5327 



5 


19683 




TABLE D.2: The ratio of the non-zero terms in the third-order cumulant function 





2. Both p and I are needed to compute the cross-correlation sequences 

Since the maximum value of d equals Nt and the maximum value of U is the 
maximum component in /, another array of the cross-correlation sequences 
can be formed such that its entry is r,(j) for i = and 

j = 0,1,2,*- - ,max{l}. The row corresoponds to the sequence of order i 
while the j*^ column corresponds to the value of the sequences at lag value 
j. Instead of using (3.96) to compute each ri{j) we use the recursion relation 
(3.92) 

n{j) = pfiij - 1 ) + P'f'i-lij) - “ 1 ) 

In this case we need to compute the values of ri(j) along the first row and 
the first column only and then compute the rest of the values recursively. The 
values in the first row are simply 

(^- 1 ) 


while the entries of the first column are given by 


ri(0) 


1 for i = 0 
0 otherwise 


(D.2) 


B. PROCEDURE FOR COMPUTING THE GRADIENT 
MATRIX H 


In the optimization algorithm to be discussed shortly, it is necessary to compute 
the gradients of the vector with respect to the model parameters. The gradient 
matrix H is thus defined as 

H = [VcX^ V^,X^ V,X„] (D.3) 

.41though this matrix can be formed from (5.21) it is preferable to avoid the unnec¬ 
essary computations by utilizing the same information stored and used to compute 



the vector Xm- Therefore each of the three partitions of H is obtained by modifying 
specific fields in the cumulant data array and using this modified array to compute 
the gradients. 

From (5.21) the gradient with respect to the vector c can be obtained by 

VcX„ = C^VcC®* (D.4) 


Then the matrix does not need to be recomputed; only the information obtained 
from the first field of the array of cumulant data is used to form another table with the 
same length. Each row of this table has N elements specifying the value of gradient 
of this term with respect to each of the model coefficients ca- If a coefficient ca is 
not included in the product of coefficients then the corresponding entry in the table 
is zero. If the coefficient is included then the corresponding value in the table is 




(D.5) 


Thus we need only to multiply the corresponding term by a value equal to the coef¬ 
ficient multiplicity in the first field then reduce the multiplicity bj’ one. 

Since a change in p changes terms in the cumulant expression that have the 

form 


rl{k)rZ{h)--rl{l,) 

the gradient of this quantity with respect to p is given by 


rZ{h)rZ{l2)--rl{l,) 




+ ^1 




^df (; j 

Notice that the derivative of each term is obtained by multiplying this term by 


T^dSh) r^{k) 

1^1 ^TTT + »^2-^77T + 

rdAh) rdjih) 




rW 

rd^Up). 


where the Ui are obtained from the original data array and the cross-correlation 
sequence values have already been computed and stored for a given range of orders 
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and time lags. Hence it remains to compute the values of r^(/) for the same range. We 
can use (3.92) to obtain a recursion for the derivatives of the correlation sequences: 

rUj) = Pr:(j - 1) + - rUi - 1) + ~ (D.6) 

P 

then precompute these values and store them in a separate array. In this array the 
entries of the first column, rj(0) are all zeros and the first row is given by 

^oU) = 

C. EXTENDED KALMAN FILTER ALGORITHM 

For the extended Kalman filter optimization procedure the model parameters 
are arranged in a vector form 

C = [c’’ pf (D.7) 

Then the recursive estimation of the value of the vector C that minimizes the quantity 
E^E is proceeds according to the structure of the Kalman Filter technique [39]. 
Assuming that at the recursion we have the value and the predection of 
is linear function of w'e can build our algorithm as follows 

Ci/i = C./i-l + A'i[Xp - 

Ct+l/t Ct/i 

Ki = 

Pi/i = Pi/i-i - P./i_iHiH'i-^HfPi/._i 

Pi+i/i = Pi/i (D.8) 

Where in this notation the matrix Hi is the gradient matrix in (D.3) computed for 
the model parameters vector Ci/i- We start the algorithm with proper initialization 
of C 3nd P. 
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